求助，程序出错

jiqixueyuan

小女子最近遇到点问题，急于求助．运用龙格库踏法求解出方程的响应：Ｘ（：，２），在此基础上对其做fft．
方程如下：dx(1)=x(3);
dx(2)=x(4);
dx(3)=-2*E*x(3)-x(1)-a*(x(1)^2+x(2)^2)*x(1)-b*(1-1/e)*(x(1)-f*x(2))+u*w^2*cos(w*t);
dx(4)=-2*E*x(4)-x(2)-a*(x(1)^2+x(2)^2)*x(2)-b*(1-1/e)*(f*x(1)+x(2))+u*w^2*sin(w*t)-G;
编写的幅值频率图的程序为：
clear
clc
clf
t0=0;tfinal=700;
t=[0:4:tfinal];
x0=[0.001 0 0.001 0];
options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
[T,X]=ode45('www',t,x0,options);
fs=200;          %采样频率
N=140001;           %采样点数
%采样时间序列s
x=X(:,2);  %生成信号
figure(1)%subplot(2,1,1);
plot(T,x);
xlabel('t/s');ylabel('位移(Y)');
title('时间历程图');
axis([0 700 -3 3]);
xf=fft(x,N);
mag=abs(xf);
df=fs/N;          %频率分辨率Hz
%绘制双边幅值谱
f=(0:N-1)*df;    %频域序列
figure(2)%subplot(2,1,2);
plot(f,abs(xf));
xlabel('Frequency(HZ)');
ylabel('Magnitude');
title('幅频谱图');
axis([0 100 0 120]);
%N=fs*t
hold off
由上面出来的幅值频率图感觉不象别人编写的，请问下，是不是我的频率幅值图的程序有问题，如果有请问下该如何修改，谢谢了，急用！

jiqixueyuan

上面的程序的运行结果图形如　fuzhitu：

而理论上应该要出现类似于　幅值图２　这样的图形．
不知道原因何在＞望大家帮忙解决，谢谢了！

jiqixueyuan

这个版块怎么就没人帮忙解决问题呢？
斑竹建议看精华，我也看了，但是还是不能解决问题，所以才发的帖子，希望能得到各位的帮忙，谢谢了！

zhangnan3509

我对你做的内容不了解，只看这两个图的比较。你做得是双边谱吧？是不是换成单边谱效果就好一点呢？

pengzk

我猜想原作者给出的谱图是系统稳态响应的谱图，而你是把系统瞬态和稳态相应的信号同时做Fourier变换了.建议你把仿真时间延长一点，用后半部分的系统相应做Fourier transform. 另外考虑一下频域分辨率的问题。