?? lyp.m
字號(hào):
#include <math.h>
#include "mex.h"
#include <stdio.h>
#include <stdlib.h>
#include <matrix.h>
double ** TwoArrayAlloc(int r,int c)
{ //兩維內(nèi)存動(dòng)態(tài)分配函數(shù)
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);
}
{
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[], // 輸出參數(shù)個(gè)數(shù)，及輸出參數(shù)數(shù)組
int nrhs, const mxArray *prhs[]) // 輸入?yún)?shù)個(gè)數(shù)，及輸入?yún)?shù)數(shù)組
{
double *pdata,*pE1;
int m1,tau1,N1,P1;
if (nrhs!=4) mexErrMsgTxt("需要4個(gè)參數(shù)!"); //檢查輸入?yún)?shù)的個(gè)數(shù)
// 取得輸入?yún)?shù)
pdata = mxGetPr(prhs[0]); // 時(shí)間序列（列向量）
tau1 = (int) *mxGetPr(prhs[1]); // 最小嵌入維數(shù)
m1 = (int) *mxGetPr(prhs[2]); // 最大嵌入維數(shù)
P1 = (int) *mxGetPr(prhs[3]); // 時(shí)間延遲
N1 = mxGetM(prhs[0]); // 序列長(zhǎng)度
// 為輸出變量分配內(nèi)存空間
plhs[0]=mxCreateDoubleMatrix(N1,1,mxREAL);
// 取得輸出參數(shù)指針
pE1 = mxGetPr(plhs[0]);
// 調(diào)用 C 運(yùn)算函數(shù) (該函數(shù)名不能和本文件名重名)
LYA_FUNCTION(pdata,tau1,m1,N1,P1,pE1);
// 求最大、最小和平均相點(diǎn)距離
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; // %最大相點(diǎn)距離
double min_d = 1.0e+100; // %最小相點(diǎn)距離
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);%相空間重構(gòu)
//%重構(gòu)相空間中相點(diǎn)的個(gè)數(shù)
Y=TwoArrayAlloc(m,M);
for (j=0;j<M;j++) //相空間重構(gòu)
{ 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)); // %平均相點(diǎn)距離
dlt_eps = (avg_d - min_d) * 0.02 ; // %若在min_eps～max_eps中找不到演化相點(diǎn)時(shí)，對(duì)max_eps的放寬幅度
min_eps = min_d + dlt_eps / 2 ; // %演化相點(diǎn)與當(dāng)前相點(diǎn)距離的最小限
max_eps = min_d + 2 * dlt_eps ; // %&&演化相點(diǎn)與當(dāng)前相點(diǎn)距離的最大限
// 從P+1～M-1個(gè)相點(diǎn)中找與第一個(gè)相點(diǎn)最近的相點(diǎn)位置(Loc_DK)及其最短距離DK
DK = 1.0e+100; // %第i個(gè)相點(diǎn)到其最近距離點(diǎn)的距離
Loc_DK = 2; // %第i個(gè)相點(diǎn)對(duì)應(yīng)的最近距離點(diǎn)的下標(biāo)
for (i=(P+1);i<=(M-1);i++) // i = (P+1):(M-1) %限制短暫分離，從點(diǎn)P+1開(kāi)始搜索
{ 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;
}
}
// 以下計(jì)算各相點(diǎn)對(duì)應(yīng)的李氏數(shù)保存到lmd()數(shù)組中
// i 為相點(diǎn)序號(hào)，從1到(M-1)，也是i-1點(diǎn)的演化點(diǎn)；Loc_DK為相點(diǎn)i-1對(duì)應(yīng)最短距離的相點(diǎn)位置，DK為其對(duì)應(yīng)的最短距離
// Loc_DK+1為L(zhǎng)oc_DK的演化點(diǎn)，DK1為i點(diǎn)到Loc_DK+1點(diǎn)的距離，稱為演化距離
// 前i個(gè)log2（DK1/DK）的累計(jì)和用于求i點(diǎn)的lambda值
// double sun_lmd,DK1,old_Loc_DK,old_DK
sum_lmd = 0 ; //% 存放前i個(gè)log2（DK1/DK）的累計(jì)和
for (i=2;i<=(M-1);i++) //i = 2 : (M-1) // % 計(jì)算演化距離
{ 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 ; //% 保存原最近位置相點(diǎn)
old_DK=DK;
// 計(jì)算前i個(gè)log2（DK1/DK）的累計(jì)和以及保存i點(diǎn)的李氏指數(shù)
if ((DK1 != 0)&&( DK != 0))
{ sum_lmd = sum_lmd + log(DK1/DK) /log(2);}
lmd[i-2] = sum_lmd/(i-1);
// 以下尋找i點(diǎn)的最短距離：要求距離在指定距離范圍內(nèi)盡量短，與DK1的角度最小
// double point_num,cos_sita,zjfwcs,dnew,DOT,CTH,cos_sita;
point_num = 0 ; //% &&在指定距離范圍內(nèi)找到的候選相點(diǎn)的個(gè)數(shù)
cos_sita = 0 ; //%&&夾角余弦的比較初值 ——要求一定是銳角
zjfwcs=0 ;//%&&增加范圍次數(shù)
while (point_num == 0)
{ //% * 搜索相點(diǎn)
for (j=1;j<=(M-1);j++)//j = 1 : (M-1)
{
if (fabs(j-i)<=(P-1)) //%&&候選點(diǎn)距當(dāng)前點(diǎn)太近，跳過(guò)！
{ continue; }
//%*計(jì)算候選點(diǎn)與當(dāng)前點(diǎn)的距離
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 )) // %&&不在距離范圍，跳過(guò)！
{continue; }
//%*計(jì)算夾角余弦及比較
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度的角，跳過(guò)！
{continue;}
if (CTH > cos_sita) // %&&新夾角小于過(guò)去已找到的相點(diǎn)的夾角，保留
{ 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) //%&&超過(guò)最大放寬次數(shù)，改找最近的點(diǎn)
{ DK = 1.0e+100;
for (ii=1;ii<=(M-1);ii++) //ii = 1 : (M-1)
{ if (fabs(i-ii)<=(P-1)) //%&&候選點(diǎn)距當(dāng)前點(diǎn)太近，跳過(guò)！
{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 ; //%&&擴(kuò)大距離范圍后重新搜索
cos_sita = 0;
}
}
}
//%取平均得到最大李雅普諾夫指數(shù)
//lambda_wolf=sum(lmd)/length(lmd);
*lambda_wolf=sum(lmd,M-2)/(M-2);
TwoArrayFree(Y);
free(lmd);
}
?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -