互相关函数的频域计算

1.时域计算

x1(n)与x2(n)的互相关定义如下x1(n)与x2(n)的互相关定义如下x_1(n)与x_2(n)的互相关定义如下

R(τ)=E[x1(m)x2(m+τ)]R(τ)=E[x1(m)x2(m+τ)]R(\tau)=E[x_1(m)x_2(m+\tau)]

离散信号的互相关由下式计算，结果中的R(n)长度为2∗N−1R(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)R(n)=\sum_{m=0}^{m=N-\left | n \right |-1}x_1(m)x_2(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τP(\omega)=\int_{-\infty }^{+\infty }R(\tau)e^{-j\omega\tau}d\tau

R(τ)=12π∫+∞−∞P(ω)ejωτdωR(τ)=12π∫−∞+∞P(ω)ejωτdωR(\tau)=\frac{1}{2\pi}\int_{-\infty }^{+\infty }P(\omega)e^{j\omega\tau}d\omega

P(ω)P(ω)P(\omega)为x1、x2x1、x2x_1、x_2的互功率谱，这一步是把互相关函数变换到了频域，互相关函数的傅里叶变化就是互谱密度，写成下式

P(ω)=∫+∞−∞∫+∞−∞x1(t)x2(t+τ)dt＊e−jωτdτP(ω)=∫−∞+∞∫−∞+∞x1(t)x2(t+τ)dt＊e−jωτdτP(\omega)=\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }x_1(t)x_2(t+\tau)dt＊e^{-j\omega\tau}d\tau

由交换积分性质和傅里叶变换的移位性质上式可简化成以下形式（参考时域卷积频域相乘推导）

P(ω)=F∗1(ω)F2(ω)P(ω)=F1∗(ω)F2(ω)P(\omega)=F_{1}^{*}(\omega)F_2(\omega)

这也是互谱密度的频域计算方法，时域互相关可以由上式做傅里叶逆变换得到

R(τ)=12π∫+∞−∞F∗1(ω)F2(ω)ejωτdωR(τ)=12π∫−∞+∞F1∗(ω)F2(ω)ejωτdωR(\tau)=\frac{1}{2\pi}\int_{-\infty }^{+\infty }F_{1}^{*}(\omega)F_2(\omega)e^{j\omega\tau}d\omega

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−12N-1时，取后2N−12N−12N-1个值，而频域计算卷积是取前部分的值，参考这里，这里取后部分是指分别取[0:+lag_max]和[-lag_max:0]的后部分，而在处理过程中使用fftshift调换了先后顺序，那么实际就相当于就取中间部分，如上面代码中的range