马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
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-8 09:26
提升卡
变色卡
显身卡
楼主
