互相关函数的频域计算

互相关函数的频域计算

1.时域计算

x1(n)x2(n) x 1 ( n ) 与 x 2 ( n ) 的 互 相 关 定 义 如 下

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

离散信号的互相关由下式计算,结果中的 R(n)2N1 R ( n ) 长 度 为 2 ∗ N − 1


R(n)=m=N|n|1m=0x1(m)x2(m+n) R ( n ) = ∑ m = 0 m = N − | n | − 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/pi

xc12 = zeros(2*N-1,1);
m = 0;
for i = -(N-1):N-1
    m = m+1;
    for t = 1:N
        if 0<(i+t)&&(i+t)<=N
            xc12(m) = xc12(m) + x2(t)*x1(t+i);
        end 
    end
end
xc12 = xc12/N;

验证可以看到自己循环计算得到的结果与matlab的xcorr结果相同

2.频域计算

由维纳-辛钦定理可知,随机信号的自相关函数和功率谱密度函数服从一对傅里叶变换的关系


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

R(τ)=12π+P(ω)ejωτdω R ( τ ) = 1 2 π ∫ − ∞ + ∞ P ( ω ) e j ω τ d ω

P(ω) P ( ω ) x1x2 x 1 、 x 2 的互功率谱,这一步是把互相关函数变换到了频域,互相关函数的傅里叶变化就是互谱密度,写成下式


P(ω)=++x1(t)x2(t+τ)dtejωτdτ P ( ω ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ x 1 ( t ) x 2 ( t + τ ) d t * e − j ω τ d τ

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


P(ω)=F1(ω)F2(ω) P ( ω ) = F 1 ∗ ( ω ) F 2 ( ω )

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


R(τ)=12π+F1(ω)F2(ω)ejωτdω R ( τ ) = 1 2 π ∫ − ∞ + ∞ F 1 ∗ ( ω ) F 2 ( ω ) e j ω τ 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以避免混叠,长度大于 2N1 2 N − 1 时,取后 2N1 2 N − 1 个值,而频域计算卷积是取前部分的值,参考这里,这里取后部分是指分别取[0:+lag_max]和[-lag_max:0]的后部分,而在处理过程中使用fftshift调换了先后顺序,那么实际就相当于就取中间部分,如上面代码中的range
  • 36
    点赞
  • 184
    收藏
    觉得还不错? 一键收藏
  • 17
    评论
评论 17
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值