互相關函數的頻域計算

1.時域計算

x1(n)與x2(n)的互相關定義如下x1(n)與x2(n)的互相關定義如下

R(τ)=E[x1(m)x2(m+τ)]R(τ)=E[x1(m)x2(m+τ)]

離散信號的互相關由下式計算，結果中的R(n)長度為2?N?1R(n)長度為2?N?1

R(n)=∑m=N?|n|?1m=0x1(m)x2(m+n)R(n)=∑m=0m=N?|n|?1x1(m)x2(m+n)

上代碼

x1 = [1,2,3,7,9,8]; x2 = [4,5,6,5,4,3]; N =length(x2); xc = xcorr(x1,x2,'biased'); [k,ind] = max(xc); an = acos((ind-N)/Fs*340/d)*180/pixc12 = zeros(2*N-1,1); m = 0; for i = -(N-1):N-1m = m+1;for t = 1:Nif 0<(i+t)&&(i+t)<=Nxc12(m) = xc12(m) + x2(t)*x1(t+i);end end end xc12 = xc12/N;

驗證可以看到自己循環計算得到的結果與matlab的xcorr結果相同

2.頻域計算

由維納-辛欽定理可知，隨機信號的自相關函數和功率譜密度函數服從一對傅里葉變換的關系

P(ω)=∫+∞?∞R(τ)e?jωτdτP(ω)=∫?∞+∞R(τ)e?jωτdτ

R(τ)=12π∫+∞?∞P(ω)ejωτdωR(τ)=12π∫?∞+∞P(ω)ejωτdω

P(ω)P(ω)為x1、x2x1、x2的互功率譜，這一步是把互相關函數變換到了頻域，互相關函數的傅里葉變化就是互譜密度，寫成下式

P(ω)=∫+∞?∞∫+∞?∞x1(t)x2(t+τ)dt＊e?jωτdτP(ω)=∫?∞+∞∫?∞+∞x1(t)x2(t+τ)dt＊e?jωτdτ

由交換積分性質和傅里葉變換的移位性質上式可簡化成以下形式（參考時域卷積頻域相乘推導）

P(ω)=F?1(ω)F2(ω)P(ω)=F1?(ω)F2(ω)

這也是互譜密度的頻域計算方法，時域互相關可以由上式做傅里葉逆變換得到

R(τ)=12π∫+∞?∞F?1(ω)F2(ω)ejωτdωR(τ)=12π∫?∞+∞F1?(ω)F2(ω)ejωτdω

matlab中xcorr函數計算相關就是在頻域計算的，這里用幾行代碼驗證下

x1 = [1,2,3,7,9,8,3,7]'; x2 = [4,5,6,5,4,3,8,2]'; N = length(x1)+length(x2)-1; NFFT = 64; range = NFFT/2+1-(N-1)/2:NFFT/2+1+(N-1)/2; xcorr(x1,x2) ifft(fft(x1,NFFT).*conj(fft(x2,NFFT))); r = fftshift(ifft(fft(x1,NFFT).*conj(fft(x2,NFFT)))); r = r(range)

關于這個計算，幾點需要注意：

• 互相關函數不是對稱的，xcorr(x1,x2) != xcorr(x2,x1)，而卷積計算是相等的，因此頻域計算要注意看誰取共軛，簡單記住哪個信號做參考就哪個信號取共軛，matlab的xcorr是第二個信號做參考
• 頻域相乘恢復到時域時得到的是[0~+lag_max，-lag_max~0]，而直接時域計算得到的就是[-lag_max~+lag_max]，因此想要與xcorr對應需要將逆變換后的數據后半部分移到前面來（fftshift），matlab 的xcorr函數內部也可以看到這個操作， % Keep only the lags we want and move negative lags before positive
  % lags.
  c = [c1(m2 - mxl + (1:mxl)); c1(1:mxl+1)];
• fft長度必須大于等于2N-1以避免混疊,長度大于2N?12N?1時，取后2N?12N?1個值，而頻域計算卷積是取前部分的值，參考這里，這里取后部分是指分別取[0:+lag_max]和[-lag_max:0]的后部分，而在處理過程中使用fftshift調換了先后順序，那么實際就相當于就取中間部分，如上面代碼中的range

總結

以上是生活随笔為你收集整理的互相关函数的频域计算的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。