声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 5157|回复: 8

[FFT] FFT变换后幅值谱及功率谱比较!

[复制链接]
发表于 2010-7-8 09:26 | 显示全部楼层 |阅读模式

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

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

x
我用傅利叶FFT变换做了这个简单的比较:共3步:
             1、由单边功率谱ss得到傅利叶幅值mag_1,结合随机相位角fai_1,得到傅利叶谱fm_1(单边),转换成傅利叶谱fm_1_1(双边),通过IFFT变换得到时程yy,如图1所示;
              2、由图2比较 由单边功率谱得到的傅利叶幅值谱mag_1和时程yy的傅利叶幅值谱mag_2,作图时均对mag_1、mag_2均*2/nfft,以得到真实的幅值。 由图3 比较功率谱ss与时程yy的功率谱pxx;      
              3、由上文得到的时程yy,进行FFT变换得到傅利叶谱fm_2,进而得到幅值谱mag_2,对时程yy进行功率谱估计得到功率谱pxx;        

     对结果比较后,有以下两个疑问:
         1、由图2,‘原幅值谱mag_1’正好为‘合成的幅值谱mag_2’的两倍,在作图时均*2/nfft,为什   
               么会有2倍之差呢?
         2、如图3对比看出,原功率谱ss和合成时程功率谱pxx相差极大,有数量级的差别,不知何故?

真心地希望大家能帮忙找出原因!

结果

结果

   Matlab程序如下: psd_wave_psd.m (3.69 KB, 下载次数: 18)

clear
Fs=100 ;                          %采样频率,HZ
Ts=20;                             %时程时间,S
%%%导出变量
nt=round(Fs*Ts+1);               %地震波数据长度
nfft=2^nextpow2(nt)              %大于并最接近nt的2的幂次方为FFT长度
df=Fs/nfft;                      %频率增量,HZ
ff=0:df:(nfft/2-1)*df;           %单边频率变化值,nt步,HZ
dt=1/Fs;                         %时间增量,S
tt=0:dt:(nt-1)*dt;              %时间变化值,nt步,S
omiga=2*pi*ff                    %圆频率变化值,nfft/2步,rad/s
%功率谱函数参数&&&&(场地、烈度、震级M)%%%%%%%%%%%
s0=128.69*0.0001; omigag=7.23; keceg=1.04; dd=0.013; omiga0=1.83;
%%%%%%%%%%%由功率谱函数计算功率谱
ss=zeros(1,nfft/2)                                    %单边功率谱
for ii=1:nfft/2
   Gomiga1(ii)=4*(keceg*omiga(ii)/omigag)^2;
   Gomiga2(ii)=(1-(omiga(ii)/omigag)^2)^2;
   Gomiga3(ii)=(1+Gomiga1(ii))/(Gomiga2(ii)+Gomiga1(ii));
   Gomiga4(ii)=s0/(1+(dd*omiga(ii))^2);
   Gomiga5(ii)=omiga(ii)^4;
   Gomiga6(ii)=(omiga0^2+omiga(ii)^2)^2;
   Gomiga(ii)=Gomiga3(ii)*Gomiga4(ii)*Gomiga5(ii)/Gomiga6(ii);
   ss(ii)=Gomiga(ii);                            %单边功率谱计算结果
end;
%%利用功率谱得到傅利叶幅值mag_1,结合随机相位角fai_1,得到傅利叶谱fm_1,通过IFFT变换得到时程yy
%功率谱转换为傅利叶幅值谱
       mag_1=sqrt(4*2*pi*df*ss)*nfft/2;      %傅利叶幅值谱,为什么乘nfft/2
      fai_1=rand(1,nfft/2)*2*pi;              %傅利叶相位角,0-2pi的随机相位角
      fm_1=mag_1.*exp(i*fai_1);               %傅利叶谱
       fm_1_1=[fm_1,fm_1(nfft/2:-1:1)];    %由单边变换至双边
      xg_1=ifft(fm_1,nfft);                 %IFFT变换
      yy=real(xg_1);                        %时程
  subplot(2,2,1);plot(yy);title('图1 合成时程yy');grid on;
%%%由上文得到的时程yy,进行FFT变换得到傅利叶谱fm_2,进而得到幅值谱mag_2和相位谱fai_2
%1.进行FFT变换
  fm_2=fft(yy,nfft)';                  %进行fft变换
  mag_2=abs(fm_2);                    %计算傅利叶幅值谱
  fai_2=angle(fm_2)
  subplot(2,2,2),plot(ff(1:nfft/2),mag_2(1:nfft/2)*2/nfft,'b'); %乘上后面的2/N得到正确幅值
        hold on; plot(ff(1:nfft/2),mag_1(1:nfft/2)*2/nfft,'r');
        title('图2 傅利叶幅值谱比较');legend('变换后幅值mag_2', '原幅值mag_1');grid on;
%%对时程yy进行功率谱估计
  pxx=(abs(fft(yy,nfft).^2)/nfft)
   subplot(2,2,3)
  plot(pxx(1:nfft/2),'b');
   %  [P,f]=psd(yy,nfft,nfft/2,hamming(nfft));
   %  plot(P);
  hold on;
  plot(ss,'-r');;
  title('图3 功率谱比较');legend('合成功率谱pxx', '原功率谱ss');grid on;

[ 本帖最后由 八戒 于 2010-7-8 09:35 编辑 ]
回复
分享到:

使用道具 举报

 楼主| 发表于 2010-7-10 23:13 | 显示全部楼层
经仔细比较发现
由:        fm_1_1=[fm_1,fm_1(nfft/2:-1:1)];    %由单边变换至双边
得到的双边傅利叶幅值谱fm_1_1,与fm_2不同。
如上例中:fm_1_1、fm_2均是[1*2048]的数组,fm_1_1是1~1024与1025~2048对称;而由yy作FFT变换得到的fm_2是关于1025点共轭,即2~1024与1026~2048是共轭的。
发表于 2010-10-8 10:36 | 显示全部楼层
学习到了   谢谢啊
发表于 2010-12-2 20:13 | 显示全部楼层
学习了。。。
发表于 2010-12-25 10:29 | 显示全部楼层
很好,学习了!谢谢!
发表于 2010-12-31 11:02 | 显示全部楼层
回复 1 # 八戒 的帖子

还是不明白结果为什么差那么大?乘nfft/2是做什么的啊?
发表于 2012-11-8 15:27 | 显示全部楼层
发表于 2012-11-8 20:32 | 显示全部楼层
发表于 2012-11-29 19:35 | 显示全部楼层
学习到了 谢谢
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-11 03:37 , Processed in 0.080630 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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