马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
clear
clc
close all
load D0603Normal29.mat %%%%%加载一个数据
x2=x(1:3000,2); %%%%%%由于数据点数太多,我取了第二个人前面3000个点,根据吕金虎书上建议的点数
SamplingFre=500; %%%%%%%采样频率为500HZ
%%%%%%%%%%%%%%%%%%%%%%中间这段是滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[b,a] = butter(5,0.75/(SamplingFre/2),'high'); %highpass Butterworth filter cutoff 0.75Hz,
y=filtfilt(b,a,x2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m_vector=2:5;
sigma=std(y);
r_vector=sigma/2*[1:4]; %%%%%%%以上4-6都是根据吕金虎的书取m和r值
N=length(y);
max_d=200; %%%%%定义最大延迟时间
S_mean_m_r_t=0;
S_mean_m_t=0;
delta_S_t=0;
S_mean_t=0;
checkMa=[];
tic
for t=1:max_d;%1:max_d;
data=disjoint(y,N,t); %%%%%将时间序列划分成t个子序列,data大小t行 * N/t列
for i=1:length(m_vector);
S_mean_m_r_t=0;
S_m_t=0;
for j=1:length(r_vector);
S_m_r_t=0;
m=m_vector(i);
r=r_vector(j);
for k=1:t;
Y1=data(k,:)'; %%%%%Y1是列向量
Y2=data(k,:); %%%%%Y2是行向量
C_I(1)=Correlation_Integral(1,Y1,r,t); %%%%求C(1,y,r,t)
X=reconstitution(Y2,m,t); %%%%对子序列相空间重构
C_I(2)=Correlation_Integral(m,X,r,t); %%%%求C(m,y,r,t)
S_m_r_t=(C_I(2)-C_I(1))+S_m_r_t; %%%%%%求S(m,r,t)
checkMa=[checkMa;C_I];
end
S_m_t(j)=(1/t)*S_m_r_t;
S_mean_m_r_t=S_m_t(j)+S_mean_m_r_t;
end
S_mean_m_t=S_mean_m_r_t+S_mean_m_t;
delta_S_m_t(i)=max(S_m_t)-min(S_m_t);
delta_S_t=delta_S_m_t(i)+delta_S_t;
end
%%%S_mean_t=S_mean_m_t+S_mean_t;
S_mean(t)=(1/16)*S_mean_m_t;
delta_S_meam(t)=(1/4)*delta_S_t;
S_cor(t)=delta_S_meam(t)+abs(S_mean(t));
end
subplot(3,1,1)
plot(delta_S_meam);
subplot(3,1,2)
plot(S_mean);
subplot(3,1,3)
plot(S_cor);
toc
这个程序我市根据吕金虎书上的步骤写的,但是好像运行不出正确的结果,dealt_S_m不知道什么原因越来越大,都达到几十几百的,我就不知道这个问题出在哪里(请大家帮忙一起看看)?本来这个程序能运行,但被我改来改去,现在出现很多NAN。 |