small_data_sets_lyapunov.c

用于小数据法求Lyapunov指数的计算的Matlab程序。

//
#include "mex.h"
#include "stdlib.h"
#include "math.h"
#define null 0
struct space
{
 struct space *next;
 double value;
};
void small_data_sets_Lyapunov(double *data,double *Length_data,double *m,double *tau,double *P,double *delt_t,double *LY_length,double *LY)
//data:时间序列数据
//Length_data:时间序列数据的长度
//tau:时间序列的延迟
//P:时间序列的平均周期
//delt_t:时间序列的样本周期
{
  int i,j,j1,j2,Count,index_j,max_i;
  int row,col,sum_former,former;
  double d_j_min,d_s,S_j_i,d_j_i;
  double *Q;
  int max_d_length;
  int P1,tau1,m1,Length_data1;
  
  struct space *d_y,*p_d_y1,*p_d_y2;
  struct space *d_length,*d_content;
  struct space *p_d_length1,*p_d_length2,*p_d_content1,*p_d_content2,*p_d_content3;
  struct space *temp;
  d_length=(struct space *)malloc(sizeof(struct space));
  d_content=(struct space *)malloc(sizeof(struct space));
  d_y=(struct space *)malloc(sizeof(struct space));
  p_d_length1=d_length;
  p_d_content1=d_content;
  p_d_y1=d_y;
  
  Length_data1=(int)(*Length_data);
  m1=(int)(*m);
  tau1=(int)(*tau);
  P1=(int)(*P);
  row=Length_data1-(m1-1)*tau1;
  col=m1;
  Q=(double *)calloc(row*col,sizeof(double));
  Count=0;
  for(i=0;i<row;i++)
    for(j=0;j<col;j++)
      Q[Count++]=data[i+j*tau1]; //相空间重构
     
  for(j1=0;j1<row;j1++)
  {
      d_j_min=10000;//寻找j1与j2的最短距离
      for(j2=0;j2<row;j2++)
      {
          d_s=0;//记录两点间的距离
          if(abs(j1-j2)>P1)
          {
            for(i=0;i<col;i++)
            {
               if ((Q[j1*col+i]-Q[j2*col+i])<0)
                  d_s=d_s-(Q[j1*col+i]-Q[j2*col+i]);
               else
                  d_s=d_s+(Q[j1*col+i]-Q[j2*col+i]);
            }
            if(d_s<d_j_min)
            {
                d_j_min=d_s;
                index_j=j2;
             }
           }   
      }
        if((row-index_j)<(row-j1))
          max_i=row-index_j;
        else
          max_i=row-j1;
        //链表操作
        p_d_length1->value=max_i;
        p_d_length2=(struct space *)malloc(sizeof(struct space)); 
        p_d_length1->next=p_d_length2;
        p_d_length1=p_d_length2;
        
        for(i=0;i<max_i;i++)
        {
           d_j_i=0;
           for(j=0;j<col;j++)
             d_j_i=d_j_i+(Q[(j1+i)*col+j]-Q[(index_j+i)*col+j])*(Q[(j1+i)*col+j]-Q[(index_j+i)*col+j]);
           //链表操作  
           p_d_content1->value=sqrt(d_j_i);
           p_d_content2=(struct space *)malloc(sizeof(struct space)); 
           p_d_content1->next=p_d_content2;
           p_d_content1=p_d_content2;
        }       
  }
    
 p_d_length1->next=null;
 p_d_content1->next=null;
 
 free(Q);
 temp=d_length;
 max_d_length=-1;
 while(temp->next!=null)
 {
 	if(max_d_length<(int)temp->value)
 		max_d_length=(int)temp->value;
 	temp=temp->next;	
 }
 for(i=0;i<max_d_length;i++)//对于每个i计算出所有的j的ln(S_j_i)的平均值
 {
	S_j_i=0;
	sum_former=0;
	Count=0;
	p_d_length1=d_length;
	p_d_length2=d_length;
	
	for(j=0;j<row;j++)
	{
		if(j==0)
			former=0;
		else
			{
				former=p_d_length1->value;
				p_d_length1=p_d_length1->next;
			}	
		sum_former=sum_former+former;
		if(i<(int)p_d_length2->value)
			{
				j1=0;
				p_d_content3=d_content;
				while(j1<(sum_former+i))
				 {
				   p_d_content3=p_d_content3->next;
				   j1=j1+1;
				  }
				if(p_d_content3->value!=0)
					{
						S_j_i=S_j_i+log(p_d_content3->value);
						Count=Count+1;
					} 	
			}
		p_d_length2=p_d_length2->next;
	}
	 p_d_y1->value=S_j_i/(Count*(*delt_t));
     p_d_y2=(struct space *)malloc(sizeof(struct space)); 
     p_d_y1->next=p_d_y2;
     p_d_y1=p_d_y2;
 }
 p_d_y1->next=null;
 p_d_y1=d_y;
 Count=0;
 while(p_d_y1->next!=null)
 {
 	LY[Count++]=p_d_y1->value;
 	p_d_y1=p_d_y1->next;
 }
 (*LY_length)=Count;
 free(d_length);
 free(d_content);
 free(d_y);
}
void mexFunction(
     int nlhs,
     mxArray *plhs[],
     int nrhs,
     const mxArray *prhs[]
)
{
	double *x1,*x2,*x3,*x4,*x5,*x6,*y1,*y2;
    
	int y1_mrows,y1_ncols;
	int y2_mrows,y2_ncols;

	//检查参量的个数
	if (nrhs!=6)
			mexErrMsgTxt("6 inputs required");
  
    x1=mxGetPr(prhs[0]);
	x2=mxGetPr(prhs[1]);
	x3=mxGetPr(prhs[2]);
	x4=mxGetPr(prhs[3]);
	x5=mxGetPr(prhs[4]);
	x6=mxGetPr(prhs[5]);
	
	  
  y1_mrows=1;
  y1_ncols=1;
  
  y2_mrows=1;
  y2_ncols=(*x2)-((*x3)-1)*(*x4);
  
      
  plhs[0]=mxCreateDoubleMatrix(y1_mrows,y1_ncols,mxREAL);
	plhs[1]=mxCreateDoubleMatrix(y2_mrows,y2_ncols,mxREAL);
      
	y1=mxGetPr(plhs[0]);
	y2=mxGetPr(plhs[1]);
  small_data_sets_Lyapunov(x1,x2,x3,x4,x5,x6,y1,y2);	   			
}