声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1184|回复: 1

[综合讨论] 讨论求解微分方程组的功率谱

[复制链接]
发表于 2009-4-3 14:13 | 显示全部楼层 |阅读模式

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

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

x
这是我用Welch法编的功率谱图,和文献上的出入很大,不知问题出在哪里,请各位大虾帮帮忙哦,呵呵呵呵

function ydot=liangfun(t,y,flag,a,b,c,d)
ydot=[a*(y(2)-y(1))+y(2)*y(3)*y(4);b*(y(1)+y(2))-y(1)*y(3)*y(4);
      -c*y(3)+6*y(2)+y(1)*y(2)*y(4);-d*y(4)+y(1)*y(2)*y(3)];

%%%%%%%%%%%%%%%%
clear;
Fs=160; %采样频率
a=50; b=4.3; c=13; d=20; w0=1.2; nfft=1024; T=2*pi/w0;
[t,y]=ode45('liangfun',[0:T/30:40],[0,1,0,0],[],a,b,c,d);
xn=y(:,1)';
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2],只计算一半的频率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx); plot_Pxx1=10*log10(Pxx1); plot_Pxx2=10*log10(Pxx2);
figure(1); plot(f,plot_Pxx); pause;
figure(2); plot(f,plot_Pxx1); pause;
figure(3); plot(f,plot_Pxx2);

[ 本帖最后由 ChaChing 于 2009-9-21 08:30 编辑 ]

功率谱图.doc

31.5 KB, 下载次数: 9

回复
分享到:

使用道具 举报

发表于 2009-9-21 00:32 | 显示全部楼层
标题有问题,没说清楚
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-13 19:55 , Processed in 0.076983 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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