从功率谱密度函数求自相关函数的问题

zhangluer

我先求方波的功率谱密度函数,再傅立叶反变换求其自相关函数,可是与xcorr函数求出的不一致（应该是三角波）。请问各位是何原因？谢谢！
N=1024;
t=0:0.01:10.23;
df=100/1024;    %频率分辨率
f=(0:N-1)*df;   
n=1:N;
w=hamming(1024);
W=diag(w);
x=square(2*pi*10*t);
X=x*W;
y=fft(X,1024);
p=y.* conj(y)
R=ifft(p);
[c,lags]=xcorr(x,'unbiased');
subplot(411);
plot(t,x)
title('原信号');
axis([0,12,-1.5,1.5]);
grid;
subplot(412);
plot(f,p)
title('功率谱');
xlabel('频率')；
grid;
subplot(413);
plot(t,R);
title('用傅立叶反变换求自相关函数');
grid;
subplot(414);
plot(lags/100,c);
title('用自带函数xcorr求自相关函数');
axis([0,5,-1,1]);
grid;

murhythm

不知道为什么不相等
不过 有偏的自相关估计才和功率谱是付氏变换对的。我看你求的自相关好像无偏的。

zhangluer

谢谢！

songzy41

可以用功率谱密度函数求自相关函数，因为功率谱密度函数和自相关函数互为傅里叶变换，但是楼主的程序中有错误，造成了计算结果的有误。
1，求功率谱密度时不能加窗函数；
2，求功率谱密度时要补零（为了能和xcorr结果比较，xcorr结果为2047个，把x补零后为2047长），楼主在没有补零下求出的相关函数是循环相关函数（circular autocorrelation function）；
3，如果直接用功率谱密度函数求出的自相关函数有偏的估算，所以在用xcorr时也要用'biased'。同时xcorr中对自相关函数归一化处理了，所以用功率谱密度函数求出的自相关函数后也作归一化处理。
把楼主的程序稍作修改如下，可看出用功率谱密度函数求出的自相关函数和自带函数求出的结果一样。
N=1024;
N2=N*2-1;
t=0:0.01:10.23;
df=100/N2;    %频率分辨率
f=(0:N2-1)*df;   
n=1:N;
w=hamming(1024);
W=diag(w);
x=square(2*pi*10*t);
%X=x*W;
y=fft(x,2047);
p=y.* conj(y);
R=real(ifft(p));
[c,lags]=xcorr(x,'biased');
subplot(411);
plot(t,x)
title('原信号');
axis([0,12,-1.5,1.5]);
grid;
subplot(412);
plot(f,p)
title('功率谱');
xlabel('频率');
grid;
subplot(413);
Rm=max(R);
plot(lags/100,fftshift(R)/Rm);
axis([-10,10,-1,1]);
title('用傅立叶反变换求自相关函数');
grid;
subplot(414);
plot(lags/100,c);
title('用自带函数xcorr求自相关函数');
axis([-10,10,-1,1]);
grid;
得图有

zhangluer

多谢你精心解答。我还几个疑问。1）为什么不能加窗2）三角波为什么衰减了，周期性的方波自相关函数应该是等幅振荡。3）R=real(ifft(p))和fftshift(R)/Rm取实部是什么意思，为什么要除以Rm？

songzy41

原帖由 zhangluer 于 2007-12-13 18:19 发表

                               
登录/注册后可看大图

多谢你精心解答。我还几个疑问。1）为什么不能加窗2）三角波为什么衰减了，周期性的方波自相关函数应该是等幅振荡。3）R=real(ifft(p))和fftshift(R)/Rm取实部是什么意思，为什么要除以Rm？

1，功率谱和自相关函数互为傅里叶变换的这一条特性，是在没有加窗的条件下。如果信号x，x加了窗w：y=x*w，则y的功率谱将和y的自相关函数互为傅里叶变换，而不是y的功率谱和x的自相关函数互为傅里叶变换；
2，周期性的方波自相关函数应该是等幅振荡，这是在周期性的方波无限长的条件下，当周期性的方波有限长时就变成三角波形式的衰减；
3，R=real(ifft(p))，在做ifft(p)时，由于计算的有限字长有可能IFFT后有虚部存在，所以取实部。又IFFT后从1025~2047是负时间延迟对应的值，所以要用fftshift(R)，把负时间延迟对应值放置在时间延迟为0的左边。
看一下xcorr对于biased和unbiased的说明：
     'biased'   - scales the raw cross-correlation by 1/M.
      'unbiased' - scales the raw correlation by 1/(M-abs(lags)).
其中的M是把数据的长度。除以Rm就相当于除M，在程序中可以改为fftshift(R)/N，得到一样的结果。如果想得到unbiased，则要增加一些运算，和xcorr一样要除以1/(N-abs(lags))。

zhangluer

原帖由 songzy41 于 2007-12-14 10:39 发表

                               
登录/注册后可看大图

1，功率谱和自相关函数互为傅里叶变换的这一条特性，是在没有加窗的条件下。如果信号x，x加了窗w：y=x*w，则y的功率谱将和y的自相关函数互为傅里叶变换，而不是y的功率谱和x的自相关函数互为傅里叶变换；
2， ...

实在太感谢你了！我现在明白了。

xurundong

w=hamming(1024);
W=diag(w);
两句可以省掉

mft0226

回答的真详细，受益颇多啊！

wanggang198684

很好，领教了

copyleft

thank you for your sharing

ab77977

回复 6 # songzy41 的帖子

plot(lags/100,c);为什么还要除以100？

csp8008

songzy41 发表于 2007-12-14 10:39

                               
登录/注册后可看大图

1，功率谱和自相关函数互为傅里叶变换的这一条特性，是在没有加窗的条件下。如果信号x，x加了窗w：y=x*w， ...

回答的真详细！
收货颇丰！谢谢！

songzy41

ab77977 发表于 2012-5-1 10:07

                               
登录/注册后可看大图

回复 6 # songzy41 的帖子

plot(lags/100,c);为什么还要除以100？

lags是以样点为单位，采样频率为100，除100，主要是把横坐标转为以时间为单位。