混沌时间序列分析源程序——感谢大家的帮助
因为毕业设计的原因，开始接触混沌时间序列分析，刚开始的确有些不好入门，对初学者来说一些程序自己编写挺困难，特别希望有高手指点或者有参考程序。我在振动论坛上得到了大家的热心帮助，作为一个“过来人”看到论坛上经常有人求助时间序列分析的相关程序，或者想参考，或者是急用，我相信大家也遇到过相似的问题——有些贴出来的程序往往是错误的或者无法运行。因为研究生阶段的课程很紧张，我已经很长时间没有来过振动论坛，以后更不会经常来，为了感谢振动论坛和大家的帮助，也为了帮助那些初学者更快的入门，我会在这个帖子里贴出我自己毕业设计时的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(,)
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 =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 =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);%求协方差矩阵
=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 =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 =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 =Kolmogorov_D2(X,Y,m_delt,tau)
%本函数用来计算关联维和Kolmogorov熵
%X:lnr满足线性区域的点
%Y:lnC满足线性区域的点
%m_delt:嵌入维的增量
%tau:时间延迟
%D2为关联维，K2为Kolmogorov熵序列
=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('周期—功率谱图')
= 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)平均
=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 =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);%计算主分量分布的标准方差
***************************************************************************************************************
时间序列分析的主要程序也就是这些了，希望大家先看一些相关的书籍，有一些基础的知识，然后参考这些程序，请大家在熟悉理解这些程序的基础上使用，因为本人为了程序的编写方便，部分程序得到的结果并不是最终的结果，需要大家寻找到线性区或这无标度区进行拟合，相信这是比较简单的的事情。
可能由于本人的疏忽，程序里有些小的错误，请大家见谅，也希望大家指正！ 非常O(∩_∩)O谢谢 真的。非常感谢，我一直都在混沌方面的程序，谢谢。 虽然我研究的和楼主相差较远，但依然对向你这样热心的人表示敬意。
向您致敬！！！
为您的努力工作和热情赞一个！！！ :victory: 感谢兄弟这么慷慨回复 楼主 yuling 的帖子
十分感谢楼主的慷慨解囊 楼主你太好了。像你这么无私的人，我很佩服。大家都应该学会并且懂得感恩 非常感谢！！！！ 感谢不胜！ 向你这样热心的人表示敬意!!!回复 楼主 yuling 的帖子
菜鸟新上路，请多指教 :handshake 感谢楼主！够慷慨！正需要这个 为什么用function =C_CMethod(data,max_d)计算延迟时间tau的时候，得到的是一堆ans值？