请教高手互相关函数和互功率谱之间的转化

songyinghua

互相关函数和互功率谱是一傅里叶变换对，那么通过对互相关函数做傅里叶变换就能得到互功率谱，或者对互功率谱做逆傅里叶变换得到相关函数。但是我使用matlab的一些指令，但是结果和预想的不同。
举个例子：对自相关函数做傅里叶变换
fs=100;N=1024;
n=0:N-1;
t=n/fs;
x=sin(10*pi*t)+cos(20*pi*t);y=cos(20*pi*t)+cos(5*pi*t);
[a,b]=xcorr(x,y);
Ft=abs(fft(a));
fa=abs(fft(a));
n1=length(fa);
n2=0:(n1-1);
fe=(n2/n1)*fs;
plot(fe(1:n1/2),fa(1:n1/2));grid on;
不知道结果是不是这样，请高手指教！
如果把互相关看成是卷积的话，那么时域卷积则频域相乘。频域里应该只剩下两信号的共有频率，其他的频率应很小。
互功率谱转化为互相关系数该怎么做呢，请指点！！

songzy41

原帖由 songyinghua 于 2010-4-11 16:28 发表
互相关函数和互功率谱是一傅里叶变换对，那么通过对互相关函数做傅里叶变换就能得到互功率谱，或者对互功率谱做逆傅里叶变换得到相关函数。但是我使用matlab的一些指令，但是结果和预想的不同。
举个例子：对自相关 ...

互相关函数和互功率谱是一傅里叶变换对，对模拟信号和对数字信号都是成立的。但在数字信号处理时，当信号为x和y，由xcorr(x,y)求出的长度差不多是x和y的二倍。LZ是从互相关函数做傅里叶变换求得互功率谱，同样可以从互功率谱做逆傅里叶变换求得互相关函数，具体编程如下：
fs=100;N=1024;
n=0:N-1;
t=n/fs;
x=sin(10*pi*t)+cos(20*pi*t);
y=cos(20*pi*t)+cos(5*pi*t);
[a,b]=xcorr(x,y);
subplot 211; plot(b,a,'k','linewidth',2); grid; hold on;
Ft=abs(fft(a));
fa=abs(fft(a));
n1=length(fa);
n2=0:(n1-1);
fe=(n2/n1)*fs;
subplot 212; plot(fe(1:n1/2),fa(1:n1/2),'k','linewidth',2);grid; hold on;
Fx=fft([x zeros(1,N)]);
Fy=fft([y zeros(1,N)]);
Fxy=Fx.*conj(Fy);
fb=abs(Fxy);
plot(fe(1:n1/2),fb(1:n1/2),'r');
legend('由xcorr->FFT求得','由FFT互谱求得')
title('互谱')
R1=real(ifft(fb));
R=ifftshift(R1);
subplot 211; plot(b,R(2:2*N),'r'); xlim([-100 100]); grid;
legend('由xcorr求得','由FFT互谱->IFFT求得')
title('互相关函数')
得图如下，不论互谱，或互相关函数都重合得较好。

ab77977

回复 2 # songzy41 的帖子

两者之间的误差很大，其中的原因是什么？具体在工程实际中应用相关函数时，其具体的算法如用C语言要如何编写？

hcharlie

有2个问题要请教
1）在有限长数组做相关时，楼上说了，需要约两倍的长度(2N-1)，请问另一半是原来数组的重复？书上讲这叫循环相关，与无穷域的相关是有区别的。
2）算例中fs=100;N=1024;不是整数周期，如换成fs=102.4;N=1024;整周期会不会好一点，这时循环相关就等于无穷相关了。