查看: 17743|回复: 110

[分形与混沌] 混沌时间序列分析源程序——感谢大家的帮助

  [复制链接]
发表于 2010-1-22 15:12 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
因为毕业设计的原因,开始接触混沌时间序列分析,刚开始的确有些不好入门,对初学者来说一些程序自己编写挺困难,特别希望有高手指点或者有参考程序。我在振动论坛上得到了大家的热心帮助,作为一个“过来人”看到论坛上经常有人求助时间序列分析的相关程序,或者想参考,或者是急用,我相信大家也遇到过相似的问题——有些贴出来的程序往往是错误的或者无法运行。
    因为研究生阶段的课程很紧张,我已经很长时间没有来过振动论坛,以后更不会经常来,为了感谢振动论坛和大家的帮助,也为了帮助那些初学者更快的入门,我会在这个帖子里贴出我自己毕业设计时的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);  %计算主分量分布的标准方差
***************************************************************************************************************

       时间序列分析的主要程序也就是这些了,希望大家先看一些相关的书籍,有一些基础的知识,然后参考这些程序,请大家在熟悉理解这些程序的基础上使用,因为本人为了程序的编写方便,部分程序得到的结果并不是最终的结果,需要大家寻找到线性区或这无标度区进行拟合,相信这是比较简单的的事情。
      可能由于本人的疏忽,程序里有些小的错误,请大家见谅,也希望大家指正!

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2010-2-15 22:01 | 显示全部楼层
非常O(∩_∩)O谢谢
发表于 2010-3-2 11:38 | 显示全部楼层
真的。非常感谢,我一直都在混沌方面的程序,谢谢。
发表于 2010-3-2 15:36 | 显示全部楼层
虽然我研究的和楼主相差较远,但依然对向你这样热心的人表示敬意。
发表于 2010-3-8 16:01 | 显示全部楼层

向您致敬!!!

为您的努力工作和热情赞一个!!!
发表于 2010-3-8 19:36 | 显示全部楼层
:victory: 感谢兄弟这么慷慨
发表于 2010-4-20 22:13 | 显示全部楼层

回复 楼主 yuling 的帖子

十分感谢楼主的慷慨解囊
发表于 2010-4-24 23:15 | 显示全部楼层
楼主你太好了。像你这么无私的人,我很佩服。大家都应该学会并且懂得感恩
发表于 2010-4-24 23:30 | 显示全部楼层
非常感谢!!!!
发表于 2010-5-30 21:53 | 显示全部楼层
感谢不胜!
发表于 2010-6-2 07:43 | 显示全部楼层
向你这样热心的人表示敬意!!!
发表于 2010-6-9 18:08 | 显示全部楼层

回复 楼主 yuling 的帖子

菜鸟新上路,请多指教
发表于 2010-6-10 11:09 | 显示全部楼层
:handshake 感谢楼主!够慷慨!正需要这个
发表于 2010-6-11 12:33 | 显示全部楼层
为什么用function [Smean,Sdeltmean,Scor,tau,tw]=C_CMethod(data,max_d)
计算延迟时间tau的时候,得到的是一堆ans值?
发表于 2010-6-17 14:02 | 显示全部楼层

回复 楼主 yuling 的帖子

感谢!请问楼主有Chaos Toolbox Ver.2.0工具箱吗,可否发给我啊,wf82426@163.com
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表