卷积, 相关与反卷积

年轻人，你不曾真正懂得数学，你只是习惯了而已。—冯 \cdot 诺伊曼

Young man, in mathematics you don't understand things. You just get used to them.—John von Neumann

卷积可用于描述过去作用对当前的影响，即卷积就是一个时空响应的叠加。

举个例子，一个地震发生了，地震波向外传播，要计算空间中任一点接收到的信号，就需要进行卷积。

即信号=源分布*点源作用的效果(格林函数)，*代表卷积运算，格林函数就是一个响应。信号就是要把这些响应在时空上叠加起来。

再如，过去女朋友生气的叠加对现在女朋友的影响也是一个卷积。所以，各位同学，这是要求积分的，总有一天要算总账的，少惹女友生气啊。

以上为个人理解，一家之言，且笑看之。

1 卷积通俗理解 (不严谨预警)

什么是卷积？通俗来讲：你在过去不同时刻惹女朋友生气的叠加，对女朋友现在坏心情的贡献就是卷积。

为什么呢？设当前时刻为t，惹女朋友生气用输入函数 \delta(t) 描述，女朋友心情用输出函数 h(t) 描述，假设女朋友现在不开心的结果是以前所有不开心的叠加（实际情况往往如此，这就是因果关系）。

注：Dirac函数\delta(t)为点源。则在 t_1 时刻，惹女朋友生气为 \delta(t-t_1) ，此时女朋友的心情为 h(t-t_1) 。为什么是 t-t_1 ,可这样理解，当时很生气（ t_1 时刻），随着时间的流逝，可能就不怎么生气了。用数学符号表示为：

\begin{equation}C\left(t-t_{1}\right)=\int_{t_1}^{t} \delta\left(\tau-t_{1}\right) h(t-\tau) d \tau \end{equation}=h(t-t_1) .

C(t) 代表当前时刻的坏心情，由于是对当前的影响，因此积分上限为t。

同理，男朋友肯定会在不同时刻 t_i 惹女朋友生气，表示为 \delta(t-t_i) ，与之对应的响应为 h(t-t_i) ,女朋友当前时刻的坏心情应是这些响应的叠加，则： C(t)=\begin{equation}\sum_{i}C(t-t_i) =\sum_{i} \int_{t_i}^{t} \delta\left复旦校花(\tau-t_{i}\right) h(t-\tau) d \tau\end{equation}

=\begin{equation}\int_{t_i }^{t}\left(\sum_{i} \delta\left(\tau-t_{i}\right)\right) h(t-\tau) d \tau\end{equation}.

令： \begin{equation}\sum_{i} \delta\left(\tau-t_{i}\right)\end{equation}=s(\tau) ,

则卷积定义为：

C(t)=\begin{equation}s(t) * h(t)=\int_{-\infty}^{\infty} s(\tau) h(t-\tau) d \tau\end{equation}=\sum_{i} h\left(t-t_{i}\right) .

关于Dirac函数，可参考：

2 卷积定义

设两函数为 f(\vec x) , g(\vec x) , \vec x 为n维向量，则 f(\vec x) 与 g(\vec x) 卷积定义为：

C(\vec x)=f(\vec x)\ast g(\vec x)=\int_{-\infty}^{+\infty}f(\vec y)g(\vec x -\vec y)d^{n}\vec y .

特别地，在一维情况下：

f(x)\ast g( x)=\int_{-\infty}^{+\infty}f( y)g(x - y)d y ,

在二维情况下：

f(x,y)\ast g(x,y)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(\alpha,\beta)g(x - \alpha,y-\beta)d \alpha d\beta ,

在三维情况下：

f(x,y,z)\ast g(x,y,z)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(\alpha,\beta,\gamma)g(x - \alpha,y-\beta,z-\gamma)d \alpha d\beta d\gamma .

3 卷积在数字信号处理中的应用

卷积是在数字信号处理中最常用，其目的是求某一输入信号经线性时不变系统（Linear time-invariant systems）作用后的输出信号。如下图所示，输入信号 x[n] 经线性系统作用后的输出信号为y(n) ，其值就等于输入信号 x[n] 与系统响应 h[n] 的卷积。

为什么是线性时不变系统？线性是一个非常优美的性质，其意味着可以叠加（任一个输入信号可以分解成多个脉冲信号，即用Dirac函数描述）；时不变意味着今天的作用和明天同样作用的作用效果相同，这实质是时间平移的不变性，说明能量守恒。

离散形式定义为(x与h进行卷积)：

y(k)=\sum_{j}^{}{x(j)h(k-j+1)}

4 卷积基本性质

1)s * h=h* s,

2) f *(s* h)=(f * s)* h,

3) h *(s+f)=h * s+h * f,

4) h * \delta=h\,

5) (s*h)'=s'*h=s*h',

一般地：

\frac{\partial ^p}{\partial x_i\partial x_j\partial x_k...}(f*g)=f*\frac{\partial ^p g}{\partial x_i\partial x_j\partial x_k...}

=\frac{\partial f}{\partial x_i}*\frac{\partial ^{p-1}g}{\partial x_j\partial x_k...}

=...

=\frac{\partial ^p f}{\partial x_i\partial x_j\partial x_k...}*g .

6)时域卷积定理： \mathcal {F}(s*h)=\mathcal {F}(s)\mathcal {F}(h).

频域卷积定理： \frac{1}{2 \pi}\mathcal {F}(s)*\mathcal {F}(h)=\mathcal {F}(s\cdot h).

时间域卷积，频率域乘积；频率域卷积，时间域乘积。

7）对线性算子 L ， L(f*g)=L(f)*g=f*L(g) ,类似于性质7.

注：\mathcal {F}(f) 代表对函数 f 进行傅里叶变化。

关于傅里叶变换，可参考如下链接：

证明：

1) s * h=\int_{-\infty}^{+\infty}s(y)h(x-y)dy

令 x-y=t \Rightarrow dy=-dt,y=x-t

RHS=\int_{+\infty}^{-\infty}s(x-t)h(t)(-dt)

=\int_{-\infty}^{+\infty}h(t)s(x-t)(dt)=h*s .

2) f *(s * h)=\int_{-\infty}^{+\infty}f(y)\left\{ \int_{-\infty}^{+\infty}s(t)h(x-y-t)dt \right\}dy

令 x-y-t=m \Rightarrow t=x-m-y,dt=-dm

RHS=\int_{-\infty}^{+\infty}f(y)\left\{ \int_{+\infty}^{\infty}s(x-m-y)h(m)(-dm) \right\}dy

=\int_{-\infty}^{+\infty}h(m)\left\{ \int_{-\infty}^{+\infty}f(y)s(x-m-y)dy \right\}dm

=h*(f*s)=(f*s)*h .

3) h *(s+f)=\int_{-\infty}^{+\infty}h(y)[s(x-y)+f(x-y)]dy

=\int_{-\infty}^{+\infty}h(y)s(x-y)+\int_{-\infty}^{+\infty}h(y)f(x-y)dy

=h*s+h*f .

4) h * \delta=\int_{-\infty}^{+\infty}h(y)\delta(x-y)dy=h(x) .

5) (s*h)'=\frac{\partial}{\partial x}\left\{ \int_{-\infty}^{+\infty}s(y)h(x-y)dy \right\}

=\int_{-\infty}^{+\infty}s(y)\frac{\partial h(x-y)}{\partial x} dy =s*h^{'}

(s*h)保研^{'}=(h*s)^{'}=h*s^{'}=s^{'}*h .

6) \mathcal {F}(s*h)=\int_{-\infty}^{+\infty}\left\{ \int_{-\infty}^{+\infty}s(y)h(x-y)dy \right\}e^{-ikx}dx

令 x-y=t\Rightarrow x=y+t,dx=dt

RHS=\int_{-\infty}^{+\infty}\left\{ \int_{-\infty}^{+\infty}s(y)h(t)dy\right\}e^{-ik(y+t)}dt

=\int_{-\infty}^{+\infty}s(y)e^{-iky}dy\int_{-\infty}^{+\infty}h(t)e^{-ikt}dt

=\mathcal {F}(s)\mathcal {F}(h) .

7）线性算子的证明与第五条性质类似卡吧，略。

5 相关(correlation)

卷积描述的是两个信号间的相互作用，相关描述的是两个信号之间的相似性，这包括自相关与互相关。用数学语言，相关定义为：

\begin{equation}R(t)=s(t)\star h(t)=\int_{-\infty}^{\infty} s(\tau-t) h(\tau) d\tau= \int_{-\infty}^{\infty} s(\tau) h(\tau+t) d\tau\end{equatio网易返现n}.

其特点为：两个函数没有反转，都是正的 \tau ；而卷积的话需要反转系统响应函数(将 \tau 变为 -\tau )。注意相关不具有交换性，即 s(t)\star h(t)\ne h(t)\star s(t) .

相关实例：

y1 =sin(x) , y2=cos(x) ,其定义域为[0, 8 \pi] ,采样间隔为dx=0.2，共计126个采样点。则y1与y2的波形及互相关系数函数如下图所示：（互相关点数：126+126-1=251）

说明：由上图可知，当cos(x)整个波形向右移动1.5时，与sin(x)吻合度最好，此时互相关系数最大；与之相对地，若cos(x)整个波形向左移动1.5时，与sin(x)吻合度最差，此时互相关系数最小；当cos(x)不移动时，其与sin(x)的互相关系数为0（对应着lags=0，corr =0），这表明此时两列波是“正交”的，这是傅里叶级数展开的基础。

因此，两个波求相关系数，（一个波不动，另一个波移动），当相关系数最大时，说明当这两列波发生相应的lags那么大的空间平移时（或时间平移，这取决于波是时间还是空间的函数），两列波最相似。做相关目的就是求两列波波形的相似程度。

相关可以看作是向量内积的推广（内积为两向量点乘，其物理意义为将一个向量投影到另一个向量）：

< \vec a ,\vec b>=\vec a \cdot \vec b=\left| a \right|\left| b\right| cos\theta .

内积越大，投影越大，两个向量间夹角越小，方向越一致，相似度越高。特别地，当内积为0时，两个向量是垂直的；只有当两个向量夹角为0时，内积最大（相关系数也最大）。因此，相关可反映两个向量空间的夹角。

当cos(x)的定义域为 [0,4\pi] 时，dx=0.2, 共计63个采样点。此时相关结果如下图所示：

当cos(x)只有原波形长度的一半时，有效相关lags减小，此时为126+63-1。最大相关系数比前图增大，这从物理上也容易理解（短波匹配长波相对于长波匹配长波更容易匹配上）。

当 y1=0.2sin(x) 时,y1与y2相关结果不变：

两列波相关系数计算(地震学常用)：

卷积与相关的转换：

a(t)\ast b(t)=\int_{-\infty}^{\infty} a(\tau) b(t-\tau) d \tau=\int_{-\infty}^{\infty} a(t-\tau) b(\tau) d \tau

=\int_{-\infty}^{\infty} a(-(\tau-t)) b(\tau) d \tau=a(-t)\star b(t) .

6 卷积实例

令 f(x)=xe^x , g(x)=H(x) ,则：

f(x)\ast g(x)=\int_{-\infty}^{+\infty} f(y)g(x-y)dy

=\int_{-\infty}^{+\infty} ye^yH(x-y)dy=\int_{-\infty}^{x} ye^ydy

=(x-1)e^x .

结果如下图所示：

解读：作用函数为 f(x) ,响应函数为 H(x) ,这意味着只要有一作用，马上就有一响应，且响应值为1，这意味着卷积结果就是作用函数的积分。若作用量为负，则卷积一定为负，且卷积负的更厉害（积分的效应）。当作用函数逐渐变正时，卷积会由负转正，卷积为0意味着， f(x) 从 -\infty 作用到当前（x=1）的积分为0.随着x的逐渐增大，卷积会越来越大，但永远不会超过作用函数，这是因为作用函数在 -\infty 到0为负。

7 用matlab实现卷积计算

再举一个例子，说明怎样用matlab中的conv函数实现两函数的卷积。

两函数为：

f(t)=sin(t) , g(t)=e^{-t} .

这两个函数在区间[0，t]卷积的解析解为：

下面代码实现：

%test conv %2022-08-03, code by Hsutyclc;clear;dt= 0.01;t = 0:dt:20;f1 = sin(t);f2 = exp(-t);f3 = 0.5*(exp(-t)+sin(t)-cos(t));fconv = conv(f2,f1)*dt;fconv_x = (1:length(fconv)).*dt+t(1);figure(1)% plot(x,f1);% hold o烟机灶具n% plot(x,f2);% hold onplot(t,f3,'-r','linewidth',2);hold onplot(fconv_x,fconv,'-.k','linewidth',1);saveas(gcf,'1.png')

在上面的代码中，注意：fconv =conv(f2,f1)*dt; 务必要乘上间隔dt，这是因为matlab中conv函数计算的是离散点的卷积，因此需要乘上时间间隔dt，才能保证振幅的正确。

红实线为解析解，黑虚线为matlab卷积结果

注意：两函数卷积后长度为：函数1点数+函数2点数-1.

两个离散点向量u，v的卷积定义为：

w(k)=\sum_{j}^{}{u(j)v(k-j+1)}

发音练习举个例子：u=[1,0,1], v=[2,7],则其卷积为：

w(1) = u(1)v(1) = 2

w(2) = u(1)v(2)+u(2)v(1) = 7

w(3) = u(1)v(3)+u(2)v(2)+u(3)v(1) =2

w(4) = u(1)v(4)+u(2)v(3)+u(3)v(2)+u(4)v(1) = 7

特点就是：翻折，移动相加(数字信号课老师可能这样讲)

离散点的卷积8 反卷积

在时间域： u(t)=f(t)*g(t)

在频率域： u(w)=f(w)g(w)

此时，若已知 u(w) 与 g(w) ,则可以利用反卷积求得 f(膝关节半月板损伤t) .

f(t)=F^{-1}\left\{ \frac{u(w)}{g(w)} \right\}

F^{-1} 代表傅里叶逆变换。即在频率域相除，再傅里叶逆变换转换到时间域。

反卷积中应特别注意的问题是：

1）g(w)中所有点的取值不能为0(不能除以0)，

2）当u(w)中含有噪声时，此时g(w)中若某些点取值很小将导致噪声影响很大，从而不能恢复出真实信号f(t),即：

u(w)=f(w)g(w)+N(w) 李丙春

\frac{u(w)}{g(w)}=f(w)+\frac{N(w)}{g(w)}

噪声影响为： \frac{N(w)}{g(w)} ，当g(w)中若某些点取值很小将导致此项很大，从而影响真实信号。

解决方法：水平线法(water-level method), PLD法(projected Landweber method)。

具体细节可参考：Bertero, M., Bindi, D., B转台轴承occacci, P., Cattaneo, M., Eva, C., & Lanza, V. (1997). Application of the projected Landweber method to the estimation of the source time function in seismology.Inverse Problems,13(2), 465.

举个例子(上面卷积的例子用来反卷积)：

8.1 不含噪声

当卷积结果不含有噪声时：

clc;clear;dt= 0.01;t = 0:dt:20;f1 = sin(t);f2 = exp(-t);f3 = 0.5*(exp(-t)+sin(t)-cos(t));fconv = conv(f2,f1)*dt;%deconvolutionf = ifft(fft(fconv)./fft([f2 zeros(1,length(f1)-1)]))/dt;%f = ifft(fft([f3 zeros(1,length(f1)-1)])./fft([f2 zeros(1,length(f1)-1)]))/dt;figure(1)plot(t,f1,'-r','linewidth',2);hold onplot(t,f(1:2001),'-k','linewidth',1);saveas(gcf,'1.png')反卷积结果非常完美，完全对上8.2 含噪声

当卷积结果含有一点点噪声时：

clc;clear;dt= 0.01;t = 0:dt:20;f1 = sin(t);f2 = exp(-t);f3 = 0.5*(exp(-t)+sin(t)-cos(t));fconv = conv(f2,f1)*dt;fconv = fconv + 0.1.*randn(size(fconv)); %add Noise%deconvolutionf = ifft(fft(fconv)./fft([f2 zeros(1,length(f1)-1)]))/dt;%f = ifft(fft([f3 zeros(1,length(f1)-1)])./fft([f2 zeros(1,length(f1)-1)]))/dt;figure(1)plot(t,f1,'-r','linewidth',2);hold onplot(t,f(1:2001),'-k','linewidth',1);saveas(gcf,'1.png')figure(2)plot(t,fconv(1:length(f3)));hold onplot(t,f3);saveas(gcf,'2.png')

此时加上噪声的卷积结果为：

蓝线为加上噪声的卷积结果，粉线为不含有噪声的真实值

此时，反卷积结果为(完全不能恢复原始信号)：

加上噪声后，反卷积失败，完全不能恢复出原始信号

为什么会这样呢？我们看看含噪声数据，不含噪声数据，以及作为分母的f2的频谱。

clc;clear;dt= 0.01;t = 0:dt:20;f1 = sin(t);f2 = exp(-t);f3 = 0.5*(exp(-t)+sin(t)-cos(t));fconv_clear = conv(f2,f1)*dt;fconv = fconv_clear + 0.1.*randn(size(fconv_clear)内蒙古生产建设兵团); %a开车打电话dd NoiseF_conv = fft(fconv);F_cle银行面试自我介绍ar = fft(fconv_clear);F_2 = fft(f2);fs = 1/dt; L =length(fconv);L1 = length(f2);f = (0:L/2)/L*fs;f1 = (0:L1/2)/L1*fs;figure(100)loglog(f,abs(F_conv(1:L/2+1)/L),'linewidth',1);hold onloglog(f,abs(F_clear(1:L/2+1)/L),'linewidth',1);hold onloglog(f1,abs(F_2(1:L1/2+1)/L1),'linewidth',1)xlabel('f(Hz)');ylabel('Spectra零售业态l amplitude');legend('Noise data','Noise-free data','denominator data');saveas(gcf,'3.png');

由此频谱可知，含噪声数据在较高频率时(约大于1 Hz)，噪声影响造成其频谱幅值远大于不含噪声时的情况，从而导致反卷积时噪声贡献非常大，因此反卷积不能恢复真实信号。

此时怎么解决呢？在该问题中噪声都是高频信号，可以先滤波再进行反卷积。而对于更一般的情况，这就需要特殊处理方法，比如水平线法，PLD方法等。PLD方法是本质上是一种梯度下降法。

含噪声的反卷积问题本质上是G矩阵为病态矩阵的线性反演问题，即：

d(w) = G(w) m(w)

已知d(w)与G(w)，求m(w)。

不过G是奇异矩阵(病态的)，该问题属于线性反演的研究范畴，在频率域利用吉洪诺夫正则化方法也能解决。

8.3 水平线法

水平线法(water level method)基本思想很简单，造成噪声影响过大的原因是分母的某些取值太小，导致相除时，噪声的影响过大。非常直观的想法就是设置一个阈值，当分母某些取值小于这个阈值时，修正此时分母的取值即可，即

（Aster, 2019, p204）

水平线法的关键是阈值的选取，需多试(结合频谱图)。

该问题用水平线法的结果为(黑线)：基本恢复。

水平线法反卷积结果(黑两性问题线)，红线为真实值c百度推广价格lc;clear;dt= 0.01;t = 0:dt:20;f1 = sin(t);f2 = exp(-t);f3 = 0.5*(exp(-t)+sin(t)-cos(t));fconv_clear = conv(f2,f1)*dt;fconv = fconv_clear + 0.1.*randn(size(fconv_clear)); %add Noisef = wldecon(fconv,f2,0.7,dt);figure(1)plot(t,f1,'-r','linewidth',2);h基金从业考试报名old onplot(t,f(1:2001),'-k','linewidth'中华传奇,1);saveas(gcf,'1.png')figure(2)plot(t,fconv(1:length(f3)));hold onplot(t,f3);saveas(gcf,'2.png')

水平线法代码为(没有经过大量测试，若有bug及时在评论区反馈)：

function [stf] = wldecon(obs,g,aerfa,dt) % water level method for deconvolution% Output:% stf: the source知识产权论文 time function% Input:% obs: seismograms of the main shock. should be a horizontal vector % g: seismograms of the aftershock, as Empirical Green's Function, or% Calculated Green's Function. also should be a horizontal vector% aerfa: it is a small number: such as 0.5,...% dt: sampling interval (s)%code by Hsuty, 2022-09-24so=size(obs);sg=size(g);if so(2)<sg(2) error('length of green function should be shorter than data!!');endF_obs = fft(obs);F_g = fft([g zeros(1,so(2)-sg(2))]);eps = max(abs(F_g))*aerfa;for ii = 1:1:length(F_o纸袋制作bs) if abs(F_g(ii)) <= eps && abs(F_g(ii)) > 0 F_g(ii) = eps*(F_g(ii)/abs(F_g(ii))); elseif abs(F_g(ii)) == 0 F_g(ii) 民权吧= eps; endendstf = ifft(F_obs./F_g)/dt;end

该问题用PLD方法，恢复的结果为(较好)：

黑线(PLD结果)，红线为真实值