lyp.m

李氏指数的计算．在ＭＡＴＬＡＢ中应用！ 程序还需要自己补充

#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[])	// 输入参数个数，及输入参数数组
{
	double *pdata,*pE1;
	int m1,tau1,N1,P1;
	
	if (nrhs!=4) mexErrMsgTxt("需要4个参数!");  //检查输入参数的个数
	
    // 取得输入参数
    pdata = mxGetPr(prhs[0]);      // 时间序列（列向量）      
    tau1 = (int) *mxGetPr(prhs[1]);            // 最小嵌入维数 
	m1 = (int) *mxGetPr(prhs[2]);            // 最大嵌入维数 
    P1 = (int) *mxGetPr(prhs[3]);      // 时间延迟   
	N1 = mxGetM(prhs[0]);              // 序列长度
	
	
    // 为输出变量分配内存空间
	plhs[0]=mxCreateDoubleMatrix(N1,1,mxREAL);
	
	
	
	// 取得输出参数指针
	pE1 = mxGetPr(plhs[0]);
    
	
    // 调用 C 运算函数 (该函数名不能和本文件名重名)
	LYA_FUNCTION(pdata,tau1,m1,N1,P1,pE1);    
	

	
 

//  求最大、最小和平均相点距离
	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);

TwoArrayFree(Y);
free(lmd);
  
    
}