|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
小女子最近遇到点问题,急于求助.运用龙格库踏法求解出方程的响应:X(:,2),在此基础上对其做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
由上面出来的幅值频率图感觉不象别人编写的,请问下,是不是我的频率幅值图的程序有问题,如果有请问下该如何修改,谢谢了,急用! |
|