声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 9388|回复: 13

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

[复制链接]
发表于 2007-12-12 22:07 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我先求方波的功率谱密度函数,再傅立叶反变换求其自相关函数,可是与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;
未命名.bmp
回复
分享到:

使用道具 举报

发表于 2007-12-12 22:26 | 显示全部楼层
不知道为什么不相等
不过 有偏的自相关估计才和功率谱是付氏变换对的。我看你求的自相关好像无偏的。

评分

1

查看全部评分

 楼主| 发表于 2007-12-13 15:31 | 显示全部楼层

回复 #2 murhythm 的帖子

谢谢!
发表于 2007-12-13 15:53 | 显示全部楼层
可以用功率谱密度函数求自相关函数,因为功率谱密度函数和自相关函数互为傅里叶变换,但是楼主的程序中有错误,造成了计算结果的有误。
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;
得图有
zhl2a.jpg

评分

1

查看全部评分

 楼主| 发表于 2007-12-13 18:19 | 显示全部楼层

回复 #4 songzy41 的帖子

多谢你精心解答。我还几个疑问。1)为什么不能加窗2)三角波为什么衰减了,周期性的方波自相关函数应该是等幅振荡。3)R=real(ifft(p))和fftshift(R)/Rm取实部是什么意思,为什么要除以Rm?
发表于 2007-12-14 10:39 | 显示全部楼层
原帖由 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))。

评分

1

查看全部评分

回复 支持 1 反对 0

使用道具 举报

 楼主| 发表于 2007-12-14 10:52 | 显示全部楼层
原帖由 songzy41 于 2007-12-14 10:39 发表

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

实在太感谢你了!我现在明白了。
发表于 2009-7-4 15:25 | 显示全部楼层
w=hamming(1024);
W=diag(w);
两句可以省掉
发表于 2011-6-8 21:32 | 显示全部楼层
回答的真详细,受益颇多啊!
发表于 2011-6-11 10:50 | 显示全部楼层
很好,领教了
发表于 2011-8-31 17:06 | 显示全部楼层
thank you for your sharing
发表于 2012-5-1 10:07 | 显示全部楼层
回复 6 # songzy41 的帖子

plot(lags/100,c);为什么还要除以100?
发表于 2013-8-22 16:09 | 显示全部楼层
songzy41 发表于 2007-12-14 10:39
1,功率谱和自相关函数互为傅里叶变换的这一条特性,是在没有加窗的条件下。如果信号x,x加了窗w:y=x*w, ...

回答的真详细!
收货颇丰!谢谢!
发表于 2013-8-22 21:07 | 显示全部楼层
ab77977 发表于 2012-5-1 10:07
回复 6 # songzy41 的帖子

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

lags是以样点为单位,采样频率为100,除100,主要是把横坐标转为以时间为单位。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 01:18 , Processed in 0.093140 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表