目录
卷积与实现
从数学的角度看,卷积就是一种运算方法,运算过程为对两个函数的计算:先翻转一个函数,再对其进行平移,然后与另一个函数相乘并积分或累加。
对于定义在连续域的函数,卷积定义为:
对于定义在离散域的函数,卷积定义为:
到此,卷积只是停留在数学计算中的一个公式。在应用中,卷积只适用于线性时不变系统 (LTI),为什么呢?(后续讨论以离散域为主)。线性就是具有可加性和齐次性:
上式中f函数为“线性系统”,a为任意常数。时不变就是当系统的输入序列存在移位时,系统响应的输出序列也有相应的移位:
因此“美好的”线性时不变系统的输出可以由它的单位脉冲响应完全描述,只要知道单位脉冲响应就能得到任意输入的输出,毕竟任何序列都可以用单位冲击序列得到。
到此,卷积有了实际应用中的“场景”——必须基于LTI系统。这个“场景”展开来看,设单位冲击序列为delta(n)的输出是单位冲击响应h(n),任意序列x(n):
因此输出y(n):
上式第一步将x(n)序列用单位冲击序列表示,第二步根据线性得到,带入delta(n)的响应,这不就是离散域的卷积定义?因此这种模样的计算公式才被定义为卷积,并且只能用于LTI系统。简记:
到此,卷积在应用中“诞生”了,因为运算实现的“场景”不同分为线性卷积、周期卷积和循环卷积。
线性卷积:两个序列乘累加计算(序列在时域上直接计算);短序列(len=M)用尾端补0的方法实现与长序列(len=N)同等长度,然后计算卷积。计算结果的长度为L=N+M-1。
周期卷积:表明周期序列与离散傅里叶系数(DFS)的关系,也就是周期序列的周期卷积对应于与之相应的离散傅里叶级数系数序列的乘积:
周期卷积与线性卷积:周期卷积是两个周期相同的序列的卷积,结果也具有相同的周期,求和是在一个周期内进行的,计算序列的周期卷积先按周期长度将序列周期延拓,短序列补0,然后在一个周期内计算卷积。线性卷积结果长度与之不同,计算过程也不同,周期卷积相当于周期移位,线性卷积是线性移位。
循环卷积:与周期卷积没有本质区别,相当于序列周期延拓后作周期卷积再取主值区间。
附上任意长度的两个复数序列直接计算卷积的python代码:
def complex_conv(A_real,A_imag,B_real,B_imag):
A_len = len(A_real)
B_len = len(B_real)
C_len = A_len + B_len - 1
C_real = [0 for i in range(C_len)]
C_imag = [0 for i in range(C_len)]
C_1 = [0 for i in range(C_len)]
C_2 = [0 for i in range(C_len)]
C_3 = [0 for i in range(C_len)]
C_4 = [0 for i in range(C_len)]
# complex conv real = conv(A_real,B_real) - conv(A_imag,B_imag)
# complex conv imag = conv(A_real,B_imag) + conv(A_imag,B_real)
for i in range(C_len):
j = max(0,i - B_len + 1)
while(j <= min(i,A_len-1)):
C_1[i] += A_real[j] * B_real[i-j]
C_2[i] += A_imag[j] * B_imag[i-j]
C_3[i] += A_real[j] * B_imag[i-j]
C_4[i] += A_imag[j] * B_real[i-j]
j += 1
C_real = [x - y for x, y in zip(C_1,C_2)]
C_imag = [x + y for x, y in zip(C_3,C_4)]
return C_real,C_imag
卷积的定义与使用条件知道后,就需要用计算去实现它。在实际的信号处理中,通常假设系统是线性时不变的,所以很多实际问题都涉及线性卷积的计算,但是根据定义式直接计算显然在实际应用中计算量过大导致耗时严重,而离散傅里叶变换的快速算法(FFT)为此提供了捷径——将线性卷积转化成循环卷积后利用傅里叶变换的性质使用FFT计算。首先傅里叶变换有时域卷积定理和频域卷积定理:
即时域的卷积对应于频域的乘积。
即时域的相乘对应于频域的卷积并除以2pi。
时域信号的线性卷积运算变换到频域的乘积运算,也就是令循环卷积的结果等于线性卷积的结果。设x(n)长度为M,h(n)长度为N,线性卷积如下:
循环卷积如下:
当循环卷积的长度L大于等于max(N,M)时,也就是信号没有被截断,并且h(n)只在区间[0,N]有非零值:
上式说明循环卷积是将线性卷积以L为周期进行周期延拓后取主值区间的结果。由于y(n)长度为(N+M-1),所以只有循环卷积的长度L大于等于线性卷积结果的长度时,这种周期延拓才不会彼此混叠,也就是循环卷积结果等于线性卷积结果,否则发生混叠二者就不能相等。此时线性卷积计算流程:将x(n)和h(n)补零延长至L,另L>=(N+M-1),DFT变换到频域相乘,IDFT反变换回时域,得到结果。而离散傅里叶变换的快速算法(FFT)大大提高了运算效率,FFT可分为按时间抽取和按频率抽取,附上MATLAB的按时间抽取的基2FFT:
%% 旋转因子生成
for i = 1:(fftnum/2)
W(i) = exp(-2i*pi*(i-1)/fftnum);
end
%% 码位倒序
j = 0;
for i = 0:(fftnum-1)
if i < j %判断是否交换数据
tmp = outdata(i+1);
outdata(i+1) = outdata(j+1);
outdata(j+1) = tmp;
tmp = dataindex(i+1);
dataindex(i+1) = dataindex(j+1);
dataindex(j+1) = tmp;
end
k = fftnum/2;
while j >= k
j = j - k;
k = k/2;
end
j = j+k;
end
%% 蝶形运算
for L = 1:M %级数
num = 2^(L-1);
for j = 0:(num-1) %各级的蝶算顺序
P = 2^(M-L)*j;
for k = j:2^L:(fftnum-2) %每次的蝶算实现
shiftres = round(outdata(k+num+1)*W(P+1)/2^15);
tmp = outdata(k+1) + shiftres;
outdata(k+num+1) = outdata(k+1) - shiftres;
outdata(k+1) = tmp;
end
end
end
IFFT可以由频域数据取共轭后进行FFT变换再乘1/N得到:
自相关与互相关
序列的相关计算:
r为互相关结果,当x和y一样时,r就是自相关结果。其定义式与卷积的定义式很像,折腾一下,相关运算就可以通过卷积运算实现。但二者的物理含义不同,线性卷积强调的是LTI系统的输入与输出的关系,相关反映的是两个信号的相似程度:
根据傅里叶变换的性质:时域的反转对应于频域的反转:
由欧拉公式得到频域的反转即取共轭,所以两个信号的相关计算可以转到频域X乘以Y的共轭,反变换回时域即可得到相关结果。