卷积、相关与FFT

设有两序列u和v，线性卷积就是把u沿y轴作镜面投影得到u(-)，然后使u(-)沿x轴正向移动，保持v不动，将两序列重叠部分求和就得到响应的卷积，如下图所示：

设size(u)=M,size(v)=N，则卷积共有M+N-1项。

循环卷积与线性卷积不同的是，u(-)和v沿圆周排列，圆周长度为L，u和v不足L的部分补0，u(-)沿顺时针转动，对重叠部分求和，显然重叠部分是整个圆周，任何时刻都是L项求和。

如果M=N=L,则循环卷积就是圆周卷积。周期卷积就是线性卷积以M+N-1为周期的周期延拓。

如果循环卷积圆周长度L=M+N-1，则线性卷积可与之等价。

时域中的线性卷积与频域中的直接相乘对应，但频域中序列的长度应为M+N-1，如果频域中长度小于M+N-1，经IFFT变换后，时域信号是原卷积信号的一段，且序列起始处有差异；循环卷积与之类似，但L的长度没有要求，所以可以利用IFFT实现循环卷积。

xn=randn(1,10);
hn=[0,4,1,2,3];
y1n=conv(xn,hn); 
Na=length(xn);
Nb=length(hn);
N=Na+Nb-1;
x1=[xn,zeros(1,N-length(xn))];
x2=[hn,zeros(1,N-length(hn))];
n=[0:N-1]; 
for i=0:N-1 
n1=mod(-n+i,N);%change serial number
y2n(i+1)=x1*x2(n1+1)';
end
y3n=ifft(fft(x1).*fft(x2));
subplot(3,1,1);stem(y1n);
subplot(3,1,2);stem(y2n); 
subplot(3,1,3);stem(y3n);

相关函数的计算算：R_xy(m)=E{x_n+my_n*}=E{x_ny_n-m*}

非归一化的计算公式为：

\( R_{xy(m)}=\sum\limits_{n=0}^{N-m-1}x_{n+m}y_n*\quad m\geqslant 0 \)

用卷积表示为：

R_xy(m)=x_(n)*y_(-n)(m)

与线性卷积类似,但相关运算中u不用翻转，只是v在u上移动，求重叠部分的和。且相关运算有性质R_xy(m)=R^*_yx(-m)，所以自相关函数是关于m=0对称的。

相应的函数为xcorr，下面为x=randn(1,100),y=xcorr(x,’biased’)的示意图：

如上图，2N-1项自相关函数R_x(n)左右对称，所以其傅里叶变换为实数，值即为功率谱密度。

由于有偏的自相关估计为：

                     \[{R_X}(n) = \frac{1}{N}x\left( { - n} \right) * x\left( n \right)\]

由于

\[x(t) = \frac{1}{{2\pi }}\int {F(\omega ){e^{i\omega t}}d\omega }  = \frac{1}{{2\pi }}\int {F(\omega )[\cos (\omega t) + i\sin (\omega t)]d\omega } \]

所以\[x( - t) = x( - t)^* = \left( {\frac{1}{{2\pi }}\int {F(\omega )[\cos (\omega t) - i\sin (\omega t)]d\omega } } \right)^* = \frac{1}{{2\pi }}\int {F^*(\omega ){e^{i\omega t}}d\omega } \]

所以

               DFT(x(-t))=DFT*(x(t))

又根据卷积与DFT的关系可得：

\[{R_X}(n) = \frac{1}{N}IFFT\left\{ {{X^*}\left( k \right)X\left( k \right)} \right\} = IFFT\left( {\frac{{{{\left| {X\left( k \right)} \right|}^2}}}{N}} \right)\]

又由于功率谱密度为自相关函数的傅里叶变换，所以功率谱密度可以表达为：

\[{P_\omega } = \frac{{{{\left| {X\left( k \right)} \right|}^2}}}{N}\]

x=cos(2*pi*5/50*[0:199])+sin(2*pi*11/50*[0:199]);  
z=fft(x,2*200-1);  
Z=z.*conj(z)/200;  
cf=ifft(Z);  
cf1=[cf(201:end),cf(1:200)];  
cx=xcorr(x,'biased');  
plot(cf1)  
hold on  
plot(cx,'r')  
hold off