|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
无水等高手们解析下,谢谢
附件是我画出的图。这是我用的程序:
close all
clear;
clc;
A=textread('d:/11136.txt');
deltt=1/512;
fs=1/deltt;
x=A(1:end,2);
t=A(1:end,1);
x=x-mean(x);
figure(1);
plot(t,x)
xlabel('Time /s');
ylabel('Displacement /mm');
title('Dispacement-time');
N=length(x);
y=zeros(1,N);
for k=1:N-1
if k==1
y(1,1)=x(1,1);
else
y(1,k)=(x(k+1,1)-x(k-1,1))./(2*deltt);
end
end
y=y';
figure(2);
%plot(x(2:5:end-1),y(2:5:end-1),'k.');grid on;
plot(x(2:1:end-1),y(2:1:end-1));grid on;
title('相图')
xlabel('x','fontsize',14);
ylabel('v','fontsize',14);
%period time
wn=[0.001 0.08];
%Filter_Coef=fir1(1000,wn,'low'); %滤波器序列,低通 ptoblem ? ?
Filter_Coef=fir1(1000,wn,'bandpass');
Filter_Coef=fir1(1000,wn,'bandpass'); %窄带阻滤波器-->陷波器
fT=filter(Filter_Coef,1,x);
%fT=conv(T,Filter_Coef);
Y=fft(fT); %进行FFT变换
Y(1)=[]; %去掉第一点,因为它是所有数据的和
N=length(Y);
power=abs(Y(1:floor(N/2))).^2/N; %求功率谱
%fs=1000/st; %st为采用时间,单位为ms
freq=(1:floor(N/2))*fs/N;
figure(3);
plot(freq,power);
title('Power spectral');
%axis([ 0 500 0 2.5e-5]);
period=1./freq;
[mp index]=max(power); %取得最大峰值的索引
max_power_freq=freq(index)
max_power_period=period(index) %得到峰值最大的周期
%庞加莱截面
%T=max_power_period; %period time
T=1/512
T0=T*2/3; % 选择截面
%phi0=pi*2/3;
for k=1:round(max(t(:,1))/T);
%for k=1:round(max(t(:,1))/2/pi);
d=t(:,1)-(k-1)*T-T0;
%d=t(:,1)-(k-1)*2*pi-phi0;
[P,K]=sort(abs(d));
x1l=x(K(1),1);
x1r=x(K(2),1);
y2l=y(K(1),1);
y2r=y(K(2),1);
t3l=t(K(1),1);
t3r=t(K(2),1);
if abs(P(1))+abs(P(2))<3e-16;
%if abs(P(1))+abs(P(2))<1;
X1(k)=x1l;
Y2(k)=y2l;
else
Q=polyfit([t3l,t3r],[x1l,x1r],1);
X1(k)=polyval(Q,(k-1)*T-T0);
%X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit([t3l,t3r],[y2l,y2r],1);
Y2(k)=polyval(Q,(k-1)*T-T0);
%Y2(k)=polyval(Q,(k-1)*2*pi-phi0);
end
end
figure(4);
%plot(X1,Y2,'k.');
plot(X1,Y2);
%plot(X1,Y2);
xlabel('x','fontsize',14);
ylabel('dx/dt','fontsize',14);
title('Poincare section'); |
|