| 因为毕业设计的原因，开始接触混沌时间序列分析，刚开始的确有些不好入门，对初学者来说一些程序自己编写挺困难，特别希望有高手指点或者有参考程序。我在振动论坛上得到了大家的热心帮助，作为一个“过来人”看到论坛上经常有人求助时间序列分析的相关程序，或者想参考，或者是急用，我相信大家也遇到过相似的问题——有些贴出来的程序往往是错误的或者无法运行。 因为研究生阶段的课程很紧张，我已经很长时间没有来过振动论坛，以后更不会经常来，为了感谢振动论坛和大家的帮助，也为了帮助那些初学者更快的入门，我会在这个帖子里贴出我自己毕业设计时的Matlab源程序以供大家参。 **************************************************************************************** function Tau=autocorrelation(data,tau_max) %data:输入的时间序列 %tau_max:最大时间延迟 %Tau:返回所求时间序列的时间延迟 N=length(data);%时间序列长度 x_mean=mean(data);%时间序列的平均值 data=data-x_mean; SSd=dot(data,data); R_xx=zeros(1,tau_max);%自相关函数初始化 for tau=1:tau_max %计算自相关函数 for ii=1:N-tau R_xx(tau)=R_xx(tau)+data(ii)*data(ii+tau); end R_xx(tau)=R_xx(tau)/SSd; end plot(1:tau_max,R_xx); hold on %作自相关函数图 line([1 tau_max],[0 0]) title('自相关法求时间延迟'); ylabel('自相关函数'); xlabel('时间延迟'); Tau=0; for jj=2:tau_max %求时间序列的时间延迟 if R_xx(jj-1)*R_xx(jj)<=0 if abs(R_xx(jj-1))>abs(R_xx(jj)) Tau=jj;break else Tau=jj-1;break end end end ****************************************************************************************************************** function C_I=correlation_integral(X,M,r) %该函数用来计算计算关联积分 %C_I:关联积分的返回值 %X:重构的相空间矢量，是一个m*M的矩阵 %M::M是重构的m维相空间中的总点数 %r:Heaviside 函数中的搜索半径 sum_H=0; for i=1:M-1 for j=i+1:M d=norm((X(:,i)-X(:,j)),inf);%计算相空间中每两点之间的距离，其中NORM(V,inf) = max(abs(V)). if r>d %sita=heaviside(r,d);%计算Heaviside 函数之值n sum_H=sum_H+1; end end end C_I=2*sum_H/(M*(M-1));%关联积分的值 ************************************************************************************************************************************************ function Data=reconstitution(data,m,tau) %该函数用来重构相空间 % m:嵌入空间维数 % tau:时间延迟 % data:输入时间序列 % Data:输出,是m*n维矩阵 N=length(data); % N为时间序列长度 M=N-(m-1)*tau; %相空间中点的个数 Data=zeros(m,M); for j=1:M for i=1:m %相空间重构 Data(i,j)=data((i-1)*tau+j); end end ***************************************************************************************** function Data=disjoint(data,t) % 此函数用于将时间序列分解成t个不相交的时间序列 % data:输入时间序列 % t:延迟，也是不相交时间序列的个数 % Data:返回分解后的t个不相交的时间序列 N=length(data); %data的长度 for i=1:t for j=1:(N/t) Data(j,i)=data(i+(j-1)*t); end end ******************************************************************************************************************************** function sita=heaviside(r,d) % 该函数用来计算Heaviside函数的值 %sita:Heaviside函数的值 %r:Heaviside函数的搜索半径 %d:两点之间的距离 if (r-d)<0 sita=0; else sita=1; end ********************************************************************************************************** function [Smean,Sdeltmean,Scor,tau,tw]=C_CMethod(data,max_d) % 本函数用于求延迟时间tau和时间窗口tw % data：输入时间序列 % max_d：最大时间延迟 % Smean，Sdeltmean,Scor为返回值 % tau：计算得到的延迟时间 % tw：时间窗口 N=length(data); %时间序列的长度 Smean=zeros(1,max_d); %初始化矩阵 Scmean=zeros(1,max_d); Scor=zeros(1,max_d); sigma=std(data); %计算序列的标准差 % 计算Smean,Sdeltmean,Scor for t=1:max_d S=zeros(4,4); Sdelt=zeros(1,4); for m=2:5 for j=1:4 r=sigma*j/2; Xdt=disjoint(data,t); % 将时间序列data分解成t个不相交的时间序列 s=0; for tau=1:t N_t=floor(N/t); % 分成的子序列长度 Y=Xdt(:,tau); % 每个子序列 %计算C(1,N/t,r,t),相当于调用Cs1(tau)=correlation_integral1(Y,r) Cs1(tau)=0; for ii=1:N_t-1 for jj=ii+1:N_t d1=abs(Y(ii)-Y(jj)); % 计算状态空间中每两点之间的距离,取无穷范数 if r>d1 Cs1(tau)=Cs1(tau)+1; end end end Cs1(tau)=2*Cs1(tau)/(N_t*(N_t-1)); Z=reconstitution(Y,m,1); % 相空间重构 M=N_t-(m-1); Cs(tau)=correlation_integral(Z,M,r); % 计算C(m,N/t,r,t) s=s+(Cs(tau)-Cs1(tau)^m); % 对t个不相关的时间序列求和 end S(m-1,j)=s/tau; end Sdelt(m-1)=max(S(m-1,:))-min(S(m-1,:)); % 差量计算 end Smean(t)=mean(mean(S)); % 计算平均值 Sdeltmean(t)=mean(Sdelt); % 计算平均值 Scor(t)=abs(Smean(t))+Sdeltmean(t); end % 寻找时间延迟tau：即Sdeltmean第一个极小值点对应的t for i=2:length(Sdeltmean)-1 if Sdeltmean(i)<Sdeltmean(i-1)&Sdeltmean(i)<Sdeltmean(i+1) tau=i; break; end end % 寻找时间窗口tw：即Scor最小值对应的t for i=1:length(Scor) if Scor(i)==min(Scor) tw=i; break; end end ***************************************************************************************************************** function [ln_r,ln_C]=G_P_good(data,tau,min_m,max_m,ss) % 本函数是利用G-P 方法计算混沌吸引子关联维 % data::待计算的时间序列 % tau: 时间延迟 % min_m:最小嵌入维 % max_m:最大嵌入维 % ss:半径搜索次数 N=length(data); %待计算的时间序列长度 ln_C=zeros((max_m-min_m)/2+1,ss); ln_r=zeros((max_m-min_m)/2+1,ss); for m=min_m:2:max_m Y=reconstitution(data,m,tau);%重构相空间 M=N-(m-1)*tau;%相空间点的数目 d=zeros(M-1,M); for i=1:M-1 for j=i+1:M d(i,j)=max(abs(Y(:,i)-Y(:,j)));%计算相空间中每两点之间的距离 end end max_d=max(max(d));%相空间中两点之间的最大距离 for i=1:M-1 %计算相空间中两点之间的最小距离 for j=1:i d(i,j)=max_d; end end min_d=min(min(d));%相空间中两点之间的最小距离 delt=(max_d-min_d)/ss;%搜索半径增加的步长 for k=1:ss r=min_d+k*delt; %C(k)=correlation_integral(Y,M,r);计算关联积分 sum_H=0; for i=1:M-1 for j=i+1:M if r>d(i,j) %计算heaviside(r,d) 函数之值 sum_H=sum_H+1; end end end C(k)=2*sum_H/(M*(M-1));%关联积分的值 ln_C((m-min_m)/2+1,k)=log(C(k)); %求lnC(r) ln_r((m-min_m)/2+1,k)=log(r); %求lnr end clear d; clear Y; plot(ln_r((m-min_m)/2+1,:),ln_C((m-min_m)/2+1,:));%画出双对数图 hold on; end **************************************************************************************************** function KL_Data=K_L_GP(data,m,tau) %该函数用来对时间序列的重构相空间进行K_L变换 %data:输入的时间序列 %m:重构相空间的维数 %tau:重构相空间的时间延迟 %KL_lamda:重构相空间协方差矩阵的m个特征值 %KL_Data:进行K_L变换后的m*m维矩阵 Data=reconstitution(data,m,tau); %对时间序列进行相空间重构 KLX=mean(Data'); %计算重构相空间矩阵每一行的平均值 KLdata=Data-KLX'*ones(1,length(data)-(m-1)*tau); %相空间矩阵每一行元素减去此行的平均值 KLData=(KLdata*KLdata')/(length(data)-(m-1)*tau); %求协方差矩阵 [V,D]=eig(KLData); %计算协方差矩阵的m个特征值和特征向量 %KL_lamda=zeros(1,m); %for ii=1:m %KL_lamda(ii)=D(ii,ii); %将协方差矩阵的特征值赋给KL_D %end KL_Data=V'*Data; %计算K_L变换后的矩阵 *********************************************************************************************************** function [ln_r,ln_C]=G_P_KL(data,tau,min_m,max_m,ss) % 本函数是利用基于KL变换的G-P 方法计算混沌吸引子关联维 % data::待计算的时间序列 % tau: 时间延迟 % min_m:最小嵌入维 % max_m:最大嵌入维 % ss:半径搜索次数 N=length(data); %待计算的时间序列长度 ln_C=zeros((max_m-min_m)/2+1,ss); ln_r=zeros((max_m-min_m)/2+1,ss); for m=min_m:2:max_m YY=K_L_GP(data,m,tau);%重构相空间并对相空间进行KL变换 Y=YY(fix(m/2.5):end,:); clear YY; M=N-(m-1)*tau;%相空间点的数目 d=zeros(M-1,M); for i=1:M-1 for j=i+1:M d(i,j)=max(abs(Y(:,i)-Y(:,j)));%计算相空间中每两点之间的距离 end end max_d=max(max(d));%相空间中两点之间的最大距离 for i=1:M-1 %计算相空间中两点之间的最小距离 for j=1:i d(i,j)=max_d; end end min_d=min(min(d));%相空间中两点之间的最小距离 delt=(max_d-min_d)/ss;%搜索半径增加的步长 for k=1:ss r=min_d+k*delt; %C(k)=correlation_integral(Y,M,r);计算关联积分 sum_H=0; for i=1:M-1 for j=i+1:M if r>d(i,j) %计算heaviside(r,d) 函数之值 sum_H=sum_H+1; end end end C(k)=2*sum_H/(M*(M-1));%关联积分的值 ln_C((m-min_m)/2+1,k)=log(C(k)); %求lnC(r) ln_r((m-min_m)/2+1,k)=log(r); %求lnr end clear d; clear Y; plot(ln_r((m-min_m)/2+1,:),ln_C((m-min_m)/2+1,:));%画出双对数图 hold on; end ********************************************************************************************************* function [ln_RS,ln_n]=Hurst(data,n_max) %本函数是用Hurst指数分析时间序列 %data：待分析的时间序列 %n_max：子序列的自大长度 %ln_RS：返回的ln(R/S)的值 %ln_n：返回的ln(n)的值 N=length(data); %待分析时间序列的长度 ln_n=log(10:10:n_max)'; %返回的ln(n)的值 for n=10:10:n_max a=floor(N/n); %时间序列分成子序列的个数 X=reshape(data(1:n*a),n,a); %把时间序列分成a个长度为n的子序列 aver=mean(X); %计算每一个子序列的平均值 cumdev=X-ones(n,1)*aver; %每个子序列的元素减去本列的平均值 cumdev=cumsum(cumdev); %计算每一个子序列的积累离差 stdev=std(X); %计算每一个子序列的均方差 RS=(max(cumdev)-min(cumdev))./stdev; %计算每一个子序列的R/S值 ln_RS(n/10,1)=log(mean(RS)); %计算所有子序列R/S值的平均值 end plot(ln_n,ln_RS) ****************************************************************************************************** function [D2,K2]=Kolmogorov_D2(X,Y,m_delt,tau) %本函数用来计算关联维和Kolmogorov熵 %X:lnr满足线性区域的点 %Y:lnC满足线性区域的点 %m_delt:嵌入维的增量 %tau:时间延迟 %D2为关联维，K2为Kolmogorov熵序列 [mm,nn]=size(X);%计算在每个嵌入维下满足线性区域的点数mm和总嵌入维数nn X_mean=mean(X);%每个嵌入维下满足线性区域的lnr平均值 Y_mean=mean(Y);%每个嵌入维下满足线性区域的lmC平均值 X_new=X-ones(mm,1)*X_mean; Y_new=Y-ones(mm,1)*Y_mean; D2=sum(sum(X_new.*Y_new))/sum(sum(X_new.*X_new));%计算关联为D2 KK=Y_mean-D2*X_mean; for ii=1:nn-1 KK_delt(ii)=KK(ii)-KK(ii+1); end K2=KK_delt/(m_delt*tau);%计算Kolmogorov熵序列 ********************************************************************************************** function T_mean=period_mean_fft(data) %该函数使用快速傅里叶变换FFT计算序列平均周期 %data：时间序列 %T_mean：返回快速傅里叶变换FFT计算出的序列平均周期 Y = fft(data); %快速FFT变换 N = length(data); %FFT变换后数据长度 Y(1) = []; %去掉Y的第一个数据，它是data所有数据的和 power = abs(Y(1:N/2)).^2; %求功率谱 nyquist = 1/2; freq = (1:N/2)/(N/2)*nyquist; %求频率 subplot(121) plot(freq,power); grid on %绘制功率谱图 xlabel('频率') ylabel('功率') title('功率谱图') period = 1./freq; %计算周期 subplot(122) plot(period,power); grid on %绘制周期－功率谱曲线 ylabel('功率') xlabel('周期') title('周期—功率谱图') [mp,index] = max(power); %求最高谱线所对应的下标 T_mean=period(index); %由下标求出平均周期 ************************************************************************************************* function lambda_1=largest_lyapunov_exponent(data,m,tau,P) 注意："这个程序得到的lambda_1不能当做最大lyapunov指数，因根据所作出的曲线选择线性区进行拟合，此处的处理是为了程序的方便" %本函数使用小数据量方法计算最大lyapunov指数 %data:时间序列 %m:嵌入维 %tau:时间延迟 %P:使用 FFT计算出的时间序列平均周期 %lambda_1:函数返回的最大最大lyapunov指数值 N=length(data); %序列长度 delt_t=1; Y=reconstitution(data,m,tau ); %重构相空间 M=N-(m-1)*tau; %M是m维重构相空间中总点数 d=zeros(M-1,M); for j=1:M d_min=1e+100; for jj=1:M %寻找相空间中每个点的最近距离点，并记下该点下标 if abs(j-jj)>P %限制短暂分离 d_s=norm(Y(:,j)-Y(:,jj));%计算分离后两点的距离 if d_s<d_min d_min=d_s; idx_j=jj; end end end max_i=min((M-j),(M-idx_j)); %计算点j的最大演化时间步长i for k=1:max_i %计算点j与其最近邻点在i个离散步后的距离 d(k,j)=norm(Y(:,j+k)-Y(:,idx_j+k)); end end %对每个演化时间步长i，求所有的j的lnd(i,j)平均 [l_i,l_j]=size(d); for i=1:l_i q=0; y_s=0; for j=1:l_j if d(i,j)~=0 q=q+1; y_s=y_s+log(d(i,j)); end end if q>0 y(i)=y_s/(q*delt_t); end end x=1:length(y); pp=polyfit(x,y,1); lambda_1=pp(1); yp=polyval(pp,x); plot(x,y,'-') ************************************************************************************************** function [Tau,I_sq]=mutual_information(data,tau_max,n) %本函数是利用互信息法求时间序列的时间延迟Tau %data：待计算的时间序列 %tau_max：最大时间延迟 %n：等间隔小格子划分数 I_sq=zeros(tau_max,1); N=length(data); %时间序列的长度 for tau=1:tau_max s=data(1:N-tau);q=data(tau+1:N); %把单变量时间序列扩充到二维相空间(s,q)上 as=min(s);bq=min(q); %在重构的相图上取框，均匀划分成n*n个小格子 delts=(max(s)-as)/n;deltq=(max(q)-bq)/n; N_sq=zeros(n); for ii=1:n %计算位于格子(ii,jj)内的相点个数 for jj=1:n for k=1:N-tau as_k=(s(k)-as)/delts; bq_k=(q(k)-bq)/deltq; if as_k>=ii-1&as_k<ii&bq_k>=jj-1&bq_k<jj N_sq(ii,jj)=N_sq(ii,jj)+1; end end end end Ntotal=sum(sum(N_sq)); Ps=sum(N_sq)/Ntotal; %计算位于一维s格子内的概率 Pq=sum(N_sq')/Ntotal; %计算位于一维q格子内的概率 Psq=N_sq/Ntotal; %计算位于二维格子(ii,jj)内概率 H_s=0; %计算s的熵 H_q=0; %计算q的熵 for i=1:n if Ps(i)~=0 H_s=H_s-Ps(i)*log(Ps(i)); elseif Pq(i)~=0 H_q=H_q-Pq(i)*log(Pq(i)); end end H_sq=0; %计算(s,q)的联合熵 for i=1:n for j=1:n if Psq(i,j)~=0 H_sq=H_sq-Psq(i,j)*log(Psq(i,j)); end end end I_sq(tau)=H_s+H_q-H_sq; %计算tau下的互信息函数 clear s q; %清空变量s和q end ********************************************************************************************** function sigma= PCA(data,m,tau) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %该函数用主分量分析(PCA)方法识别混沌和噪声。混沌信号的主分量谱图应是一条过定 %点且斜率为负值的直线，而噪声信号的主分量谱图应是一条与X轴接近平行的直线，故 %可以用主分量分布标准方差作为识别混沌和噪声的一种特征。 %data：输入的待分析时间序列 %m：重构相空间的维数 %tau:重构相空间的时间延迟 %sigma：主分量分布的标准方差 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Data=reconstitution(data,m,tau); %对时间序列进行相空间重构 KLX=mean(Data'); %计算重构相空间矩阵每一行的平均值 KLdata=Data-KLX'*ones(1,length(data)-(m-1)*tau); %相空间矩阵每一行元素减去此行的平均值 A=(KLdata*KLdata')/(length(data)-(m-1)*tau); %求协方差矩阵 %A=(Data*Data'); lamda=eig(A); %计算协方差矩阵的特征值 lamda=-sort(-lamda); %将特征值从大到小排序 lamda_PCA=log(lamda/sum(lamda)); plot(lamda_PCA); %做主分量谱图 ylabel('PCA') title('主分量谱图') sigma=std(lamda_PCA); %计算主分量分布的标准方差 *************************************************************************************************************** 时间序列分析的主要程序也就是这些了，希望大家先看一些相关的书籍，有一些基础的知识，然后参考这些程序，请大家在熟悉理解这些程序的基础上使用，因为本人为了程序的编写方便，部分程序得到的结果并不是最终的结果，需要大家寻找到线性区或这无标度区进行拟合，相信这是比较简单的的事情。 可能由于本人的疏忽，程序里有些小的错误，请大家见谅，也希望大家指正！ |
GMT+8, 2024-11-26 12:01 , Processed in 0.070269 second(s), 23 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.
