打开APP
userphoto
未登录

开通VIP,畅享免费电子书等14项超值服

开通VIP
互相关函数的频域计算

互相关函数的频域计算

1.时域计算

x1(n)x2(n)

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

离散信号的互相关由下式计算,结果中的R(n)2N1


R(n)=m=N|n|1m=0x1(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-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 endendxc12 = xc12/N;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18

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

2.频域计算

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


P(ω)= R(τ)ejωτdτ

R(τ)=12π P(ω)ejωτdω

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


P(ω)= x1(t)x2(t τ)dtejωτdτ

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


P(ω)=F1(ω)F2(ω)

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


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)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9

关于这个计算,几点需要注意:

  • 互相关函数不是对称的,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时,取后2N1个值,而频域计算卷积是取前部分的值,参考这里,这里取后部分是指分别取[0: lag_max]和[-lag_max:0]的后部分,而在处理过程中使用fftshift调换了先后顺序,那么实际就相当于就取中间部分,如上面代码中的range
本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报
打开APP,阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
用MATLAB对WAV文件做FFT、IFFT及短时傅里叶变换code
傅里叶变换的Matlab代码与注释
炒股公式FFT
SSVEP脑机接口及数据集处理
MATLAB与信号处理_3
画包络谱和幅值谱matlab函数示例代码
更多类似文章 >>
生活服务
热点新闻
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权!
如果VIP功能使用有故障,
可点击这里联系客服!

联系客服