八戒 发表于 2010-7-8 09:26

FFT变换后幅值谱及功率谱比较！

我用傅利叶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程序如下：

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=;    %由单边变换至双边
      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');
   %=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_1，与fm_2不同。
如上例中：fm_1_1、fm_2均是的数组，fm_1_1是1~1024与1025~2048对称；而由yy作FFT变换得到的fm_2是关于1025点共轭，即2~1024与1026~2048是共轭的。

lidaijin 发表于 2010-10-8 10:36

学习到了   谢谢啊

robustcococole 发表于 2010-12-2 20:13

学习了。。。

sunlm01 发表于 2010-12-25 10:29

很好，学习了！谢谢！

glwh 发表于 2010-12-31 11:02

回复 1 # 八戒 的帖子

还是不明白结果为什么差那么大？乘nfft/2是做什么的啊？

电光幻影 发表于 2012-11-8 15:27

{:{39}:}

kyu16866 发表于 2012-11-8 20:32

glwh 发表于 2010-12-31 11:02 static/image/common/back.gif
回复 1 # 八戒 的帖子

还是不明白结果为什么差那么大？乘nfft/2是做什么的啊？

http://forum.chinavib.com/thread-49413-1-1.html
参考《看看这里有没有你要问的问题——信号处理专栏话题索引》

午餐的雨 发表于 2012-11-29 19:35

学习到了 谢谢
页: [1]
查看完整版本: FFT变换后幅值谱及功率谱比较！