卷積、相關與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