互相关函数的时域和频域计算(matlab)

您所在的位置:网站首页 matlab互相关函数并画图 互相关函数的时域和频域计算(matlab)

互相关函数的时域和频域计算(matlab)

2024-07-08 19:41| 来源: 网络整理| 查看: 265

       相关是平时工作用到很多的概念,平时老是不经意会犯一些细节上的错,比如得到的顺序,频域做出来不一致等,所以写这篇文章吧这些细节好好扣了一下。

相关的理解

        说两个东西相关,可以说两个东西的相似程度, 在数学上抽象化这一特征得出来的公式是这样的:

                                                   r_{xy}[k] = \sum_{n} x[n] \cdot y[n+k]

        其中r_{xy}[k] 代表当 y[n]相对x[n]滑动 k个单位时的互相关值。

        为什么这个式子可以表征两个东西的相似程度呢?我是这么理解的:每一个元素都代表对象的一个属性,如果两个对象的同一个属性为正数或者都为负数(两个对象在这个属性上类似)那么乘起来加到累加和里面是会增大这个相关值的,相反则会减小这个相关值。

         当自己与自己做相关时(不滑动)可以获得最大值,这里拿一个最简单的序列a,b来说,a,b与a,b的相关值肯定会大于a,b与b,a的相关值,因为a^2+b^2>2ab。扩展为更长的序列也是这样,永远是自己与自己做相乘累加可以获得最大值。

        在实际应用中,比如雷达检测目标的过程中,波形发送出去再回来后,它本身会发生时延,以及幅值也会随着噪声发生变化,但我们通过用原来的信号去和这个接收到的信号做相关得到相关值函数,函数中一般能获得一个最大值,通过这个最大值所在的位置就可以获得目标的位置。

时域计算

        在时域中计算相关,matlab提供了xcorr函数,它实际上就是把一个序列固定A,另一个序列B从最后一位对齐序列A的第一位到序列B的第一位对齐序列A的最后一位,每一次移动一位的同时再将对应的值相乘再累加。

代码如下:

A = [1, 2, 3]; % 信号A B = [4, 5, 6]; % 信号B lenA = length(A); lenB = length(B); R = zeros(1, lenA+lenB-1); % 在A上添加0,使其长度与结果相同 A_add0 = [zeros(1, lenA-1), A, zeros(1, lenA-1)]; % 计算相关 for i = 1:(lenA + lenB - 1) % A与B当前对齐的部分相乘后求和 R(i) = sum(B .* A_add0(i:i+lenA-1)); end

        大家可以自行验证一下和xcorr计算的一致性。这里有一点要注意一下,xcorr(A,B)输出的第一个数是B的最右边和A的最左边相乘的结果,最后一个数是B的最左边和A的最右边相乘的结果。

频域计算

        频域的相乘等于时域的卷积,时域的卷积和相关不同的是,它计算时需要把序列反转再去做相乘累加。

        傅立叶变换在处理信号时具有一个重要的性质:对信号取共轭复数在时间域相当于时间反转(即f(t)变为f(-t))。那么只要我们做频域相乘的时候把其中一个取共轭,就可以得到时域的相关。

        这里还涉及到一个循环卷积和线性卷积的问题:(我们做的是相关,但为了方便我直接说卷积了,相关是一样的)直接把两个信号做FFT,取共轭相乘,再做IFFT得出来的是循环卷积的结果。

        刚刚我们在时域做相关的时候,第一个数是B的最右边和A的最左边相乘的结果,也就是这样

                                             4       5       6

                                                               1      2     3

        但用FFT等效的循环卷积,它会把超出相乘范围的值移动到另一边去,当算到B的最右边和A的最左边相乘时,本来其他位置应该用0计算的,却成了序列中其他的元素去计算了:

                                                               6      4     5

                                                               1      2     3

        所以为了让我们的结果能正确,我们需要给两个序列补0,使得每个序列长度为len(A)+len(B)-1,差多少补多少。这里A和B的长度都为3,所以要把他们都补到3+3-1=5;此时

 A=[1 2 3 0 0 ],B=[4 5 6 0 0],我们再来看这个时候的相乘累加值                               

                                                 6      0      0      4     5     

                                                 1      2      3      0     0

        可以看到,我们通过补0来使得循环卷积获得线性卷积一样的结果。

        代码如下:

N = length(A)+length(B)-1; % 计算A和B的FFT FA = fft(A, N); FB = fft(B, N); % 乘以B的共轭复数FFT并做IFFT得到互相关,fftshift完成频谱搬移 r0= fftshift(ifft( FA .*conj(FB), N));

当然,用频域计算得到的结果的顺序和时域不太一样。

横坐标顺序问题

        对于信号A(长度为 M)和信号B(长度为 N),不指定额外的参数时,xcrorr(A,B)函数返回一个长度为2*max(M,N)-1的向量,其中包含了所有可能的滞后值的互相关。例如,如果M=16和N=32,则滞后范围从-31​ 到+31。

        关于频域计算的结果,通过将输入信号的快速傅立叶变换(FFT)相乘来计算互相关,然后执行逆 FFT(IFFT)来回到时间域。在没有进行任何移位之前,你会得到一个以0开始的序列,时间滞后从0增加到正方向,表示信号之间的相关性。

        要令频域方法求得的互相关结果与xcrorr相对应,需要调整频域结果的顺序。这可以通过循环移位完成,将结果向左移动 ​2*max(M,N)-1个元素,使得零滞后点位于输出数组的中央。进行该移位后,你得到的频域互相关结果将具有与 ​xcrorr 相同的滞后值排序,由负滞后到正滞后排列。

实例

        分别用时域和频域的计算结果如图:

       实例代码如下:

% 假设 x 和 y 是你的两个输入信号 x =[ 0 45 63 45 0 -45 -63 -45 0 -45 -63 -45 0 45 63 45 ]; % 举例的信号 x[n] y =[ 9 6 14 0 7 6 7 12 5 12 52 64 48 11 -38 -61 -40 9 -42 -52 -41 14 49 74 48 4 1 9 10 8 6 10 ]; % % 举例的信号 y[n] % 使用 xcorr 获取时间域的互相关 [R,lags] = xcorr(x, y); % 解释 xcorr 的结果排序 % lags = -(numel(y)-1):(numel(x)-1); % 这是 xcorr 生成的滞后值序列 % 进行频域计算 X = fft(x, 63); Y = fft(y, 63); R_freq = ifft(X .* conj(Y)); % 循环移位以使 0 滞后对齐于输出向量的中心位置 R_freq_shifted = circshift(R_freq, [0, numel(y)-1]); % 补充说明,循环移位后,R_freq_shifted 的排序现在类似于 xcorr 输出的排序 figure,plot(lags,R,'r'),hold plot(lags,R_freq_shifted,'b.') xlabel('时延') ylabel('相关值幅度') legend('xcorr计算结果0','频域计算结果')

        注意这里FFT计算点数为63(2max(N,M)-1),主要是为了和xcorr保持一致。上一节讲到FFT点数(补零后的点数)至少为N+M-1是为了不让循环卷积产生混叠,但要覆盖所有的移位情况的话还是得取2max(N,M)-1。

写在最后

        没有把所有细节都写上去,比如共轭取哪个,相乘的顺序等也会有影响。不去深究的话,按照文章中代码顺序也可以得到正确的结果了,希望对大家有所帮助~



【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3