?? lyap_wolf.cpp
字號:
#include <math.h>
//#include "mex.h"
#include <stdio.h>
#include <stdlib.h>
//#include <matrix.h>
double ** TwoArrayAlloc(int r,int c)
{ //兩維內存動態分配函數
double *x, **y;
int n;
x=(double *)calloc(r*c,sizeof(double));
if(x==NULL)
{
printf("\nFAilue in memory applying.");
//AfxMessageBox("Failure in memory applying.");
exit(0);
}
y=(double **)calloc(r,sizeof(double *));
for(n=0;n<=r-1;++n)
y[n]=&x[c*n];
return(y);
}
void TwoArrayFree(double **x)//2 維內存釋放函數
{
free(x[0]);
free(x);
}
double sum(double *array,int n)
{
int i;
double result=0.0;
for (i=0;i<n;i++)
result=result+array[i];
return result;
}
/*
void mexFunction (int nlhs, mxArray *plhs[], // 輸出參數個數，及輸出參數數組
int nrhs, const mxArray *prhs[]) // 輸入參數個數，及輸入參數數組
{
int i;
if (nrhs!=1) mexErrMsgTxt("只需要1個參數!"); //檢查輸入參數的個數
// 取得輸入參數
pdata = mxGetPr(X); // 時間序列（列向量）
length = mxGetM(X); // 序列長度
// for(i=0;i<length;i++)
//printf("%d %f\n",i ,pdata[i]);
// 為輸出變量分配內存空間
//T= mxCreateDoubleMatrix(1,1,mxREAL); //用于存放E1
E= mxCreateDoubleMatrix(101,1,mxREAL); //用于存放E2
// T= (int*)malloc(sizeof(int)); //用于存放E1
//E=(double*)malloc((corrlength+1)*sizeof(double)); //用于存放E2
// 取得輸出參數指針
//tau = mxGetPr(T);
entropy = mxGetPr(E);
// 調用 C 運算函數 (該函數名不能和本文件名重名)
//tau=(int*)malloc(sizeof(int));
// entropy=(double*)malloc((corrlength+1)*sizeof(double));
mutual_FUNCTION(pdata,length,partitions,corrlength,array,h1,h11,h2);
//free(pdata);free(entropy);free(array);free(h1);free(h11);free(h2);
// for(i=0;i<corrlength;i++)
//printf("%d %f\n",i ,entropy[i]);
//printf("tau=%d\n",*tau);
}
*/
void main()
{
int qq,length=0;
//double *series,min,interval,shannon;
FILE *file,*fp;
char a[25];
double *pdata,ndata;//result,*value;
if(!(fp=fopen("test.txt","r")))
{
printf("打開文件數據錯誤!\n");
exit(0);
}
//得到數據個數 size
while(fscanf(fp,"%f",&ndata)==1)
{length++;}
rewind(fp);
//Set pointer to beginning of file:
fseek( fp, 0L, SEEK_SET );
//初始化數據
pdata=(double*)malloc(length*sizeof(double));
//Read data back from file:
for(qq=0;qq<length;qq++)
{ fgets(a,25,fp);
pdata[qq]=atof(a);}
fclose( fp );
/*
file=fopen("file_out.txt","w");
for (k=0;k<=corrlength;k++)
{ fprintf(file,"%d %e\n",k,entropy[k]);}
fprintf(file,"tau=%d",*tau);
fclose(file);*/
// function lambda_wolf=lyapunov_wolf(data,N,m,tau,P)
/*
% 該函數用來計算時間序列的最大Lyapunov 指數--Wolf 方法
% m: 嵌入維數
% tau:時間延遲
% data:時間序列
% N:時間序列長度
% P:時間序列的平均周期,選擇演化相點距當前點的位置差，即若當前相點為I，則演化相點只能在|I－J|>P的相點中搜尋
% lambda_1:返回最大lyapunov指數值
min_point=1 ; %&&要求最少搜索到的點數
MAX_CISHU=5 ; %&&最大增加搜索范圍次數
%FLYINGHAWK*/
// 求最大、最小和平均相點距離
int i,j,k,ii,tau=6,m=8,P=50,M,N,Loc_DK,old_Loc_DK;
double **Y,*lmd;
double max_d = 0,d; // %最大相點距離
double min_d = 1.0e+100; // %最小相點距離
double avg_dd = 0;
double dlt_eps,min_eps,max_eps,DK;
double sum_lmd,DK1,old_DK,avg_d;
double point_num,cos_sita,zjfwcs,dnew,DOT,CTH,*lambda_wolf;
int min_point=1;
int MAX_CISHU=5;
N=length;
M=N-(m-1)*tau;
lambda_wolf=(double*)malloc(sizeof(double));
lmd=(double*)malloc((M-2)*sizeof(double));
//Y=reconstitution(data,N,m,tau);%相空間重構
//%重構相空間中相點的個數
Y=TwoArrayAlloc(m,M);
for (j=0;j<M;j++) //相空間重構
{ for (i=0;i<m;i++)
Y[i][j]=*(pdata+i*tau+j);
}
for (i=1;i<=(M-1);i++) //i = 1 : (M-1)
{ for (j=i+1;j<=M;j++) //j = i+1 : M
{ d = 0;
for (k=1;k<=m;k++) //k = 1 : m
{ d = d + (Y[k-1][i-1]-Y[k-1][j-1])*(Y[k-1][i-1]-Y[k-1][j-1]); }
d = sqrt(d);
if (max_d < d)
{ max_d = d;}
if (min_d > d)
{ min_d = d;}
avg_dd = avg_dd + d;
}
}
avg_d = 2*avg_dd/(M*(M-1)); // %平均相點距離
dlt_eps = (avg_d - min_d) * 0.02 ; // %若在min_eps～max_eps中找不到演化相點時，對max_eps的放寬幅度
min_eps = min_d + dlt_eps / 2 ; // %演化相點與當前相點距離的最小限
max_eps = min_d + 2 * dlt_eps ; // %&&演化相點與當前相點距離的最大限
// 從P+1～M-1個相點中找與第一個相點最近的相點位置(Loc_DK)及其最短距離DK
DK = 1.0e+100; // %第i個相點到其最近距離點的距離
Loc_DK = 2; // %第i個相點對應的最近距離點的下標
for (i=(P+1);i<=(M-1);i++) // i = (P+1):(M-1) %限制短暫分離，從點P+1開始搜索
{ d = 0;
for (k=1;k<=m;k++) //k = 1 : m
{ d = d + (Y[k-1][i-1]-Y[k-1][0])*(Y[k-1][i-1]-Y[k-1][0]);}
d = sqrt(d);
if ((d<DK)&&(d>min_eps))
{ DK = d;
Loc_DK = i; //Loc_DK = i;
}
}
// 以下計算各相點對應的李氏數保存到lmd()數組中
// i 為相點序號，從1到(M-1)，也是i-1點的演化點；Loc_DK為相點i-1對應最短距離的相點位置，DK為其對應的最短距離
// Loc_DK+1為Loc_DK的演化點，DK1為i點到Loc_DK+1點的距離，稱為演化距離
// 前i個log2（DK1/DK）的累計和用于求i點的lambda值
// double sun_lmd,DK1,old_Loc_DK,old_DK
sum_lmd = 0 ; //% 存放前i個log2（DK1/DK）的累計和
for (i=2;i<=(M-1);i++) //i = 2 : (M-1) // % 計算演化距離
{ DK1 = 0.0;
for (k=1;k<=m;k++) //k = 1 : m
{ // DK1 = DK1 + (Y[k][i]-Y[k][Loc_DK+1])*(Y[k][i]-Y[k][Loc_DK+1]);
DK1= DK1+(Y[k-1][i-1]-Y[k-1][Loc_DK])*(Y[k-1][i-1]-Y[k-1][Loc_DK]);
//DK1= DK1+(*(pdata+(k-1)*tau+i-1)-*(pdata+(k-1)*tau+i))*(*(pdata+(k-1)*tau+i-1)-*(pdata+(k-1)*tau+i));
}
DK1 = sqrt(DK1);
old_Loc_DK = Loc_DK ; //% 保存原最近位置相點
old_DK=DK;
// 計算前i個log2（DK1/DK）的累計和以及保存i點的李氏指數
if ((DK1 != 0)&&( DK != 0))
{ sum_lmd = sum_lmd + log(DK1/DK) /log(2);}
lmd[i-2] = sum_lmd/(i-1);
// 以下尋找i點的最短距離：要求距離在指定距離范圍內盡量短，與DK1的角度最小
// double point_num,cos_sita,zjfwcs,dnew,DOT,CTH,cos_sita;
point_num = 0 ; //% &&在指定距離范圍內找到的候選相點的個數
cos_sita = 0 ; //%&&夾角余弦的比較初值 ——要求一定是銳角
zjfwcs=0 ;//%&&增加范圍次數
while (point_num == 0)
{ //% * 搜索相點
for (j=1;j<=(M-1);j++)//j = 1 : (M-1)
{
if (fabs(j-i)<=(P-1)) //%&&候選點距當前點太近，跳過！
{ continue; }
//%*計算候選點與當前點的距離
dnew = 0;
for (k=1;k<=m;k++)//k = 1 : m
{ dnew = dnew + (Y[k-1][i-1]-Y[k-1][j-1])*(Y[k-1][i-1]-Y[k-1][j-1]);}
dnew = sqrt(dnew);
if ((dnew < min_eps)||( dnew > max_eps )) // %&&不在距離范圍，跳過！
{continue; }
//%*計算夾角余弦及比較
DOT = 0;
for (k=1;k<=m;k++)//k = 1 : m
{ DOT = DOT+(Y[k-1][i-1]-Y[k-1][j-1])*(Y[k-1][i-1]-Y[k-1][old_Loc_DK]);}
CTH = DOT/(dnew*DK1);
if (acos(CTH) > (3.14151926/4) ) // %&&不是小于45度的角，跳過！
{continue;}
if (CTH > cos_sita) // %&&新夾角小于過去已找到的相點的夾角，保留
{ cos_sita = CTH;
Loc_DK = j;
DK = dnew;
}
point_num = point_num +1;
}
if (point_num <= min_point)
{
max_eps = max_eps + dlt_eps;
zjfwcs =zjfwcs +1;
if (zjfwcs > MAX_CISHU) //%&&超過最大放寬次數，改找最近的點
{ DK = 1.0e+100;
for (ii=1;ii<=(M-1);ii++) //ii = 1 : (M-1)
{ if (fabs(i-ii)<=(P-1)) //%&&候選點距當前點太近，跳過！
{continue;}
d = 0;
for (k=1;k<=m;k++)//k = 1 : m
{ d = d + (Y[k-1][i-1]-Y[k-1][ii-1])*(Y[k-1][i-1]-Y[k-1][ii-1]);}
d = sqrt(d);
if ((d<DK)&&(d>min_eps))
{ DK = d;
Loc_DK = ii;
}
}
break;
}
point_num = 0 ; //%&&擴大距離范圍后重新搜索
cos_sita = 0;
}
}
}
//%取平均得到最大李雅普諾夫指數
//lambda_wolf=sum(lmd)/length(lmd);
*lambda_wolf=sum(lmd,M-2)/(M-2);
printf("lambda1=%f\n",*lambda_wolf);
TwoArrayFree(Y);
free(lmd);
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -