das.c

用小数据量方法求lyapunov指数

//
#include "mex.h"
#include "stdlib.h"

double heaviside(double x,double y)
{
  if(x>=y)
   return 1.0;
  else
   return 0.0;
}
double correlation_integ(int m,int t,double r,int size_col,double *Q_temp)
{
 int i,j,k,L;
 int row,col;
 double *Q,S,max_S,Sum;
 
 row=size_col-(m-1)*t;
 if(row<=1)
   return 0.0;
 col=m;
 k=0;
 Sum=0;
 Q=(double *)calloc(row*col,sizeof(double));
 for(i=0;i<row;i++)
     for(j=0;j<col;j++)
       Q[k++]=Q_temp[i+j*t]; //把Q_temp进行重构相空间，存到Q中
 for(i=0;i<row;i++)
    {
      for(j=i+1;j<row;j++)
        {
          max_S=-100;
          for(L=0;L<col;L++)
         {
           S=Q[i*col+L]-Q[j*col+L];
           //将S转为|S|
           if(S<0)
             S=-S;
           if(S>max_S)
             max_S=S;
         }
         Sum=Sum+heaviside(r,max_S);
       }
    } 
  Sum=2*Sum/(row*(row-1));    
 free(Q);
 return Sum;
 }


void DAS(double *L_Data,double *Data,double *t_max,double *R,double *delt_ave_s,double *ave_s,double *S_cor)
{
  int t,i,j,i1,j1;
  double *Q,*Q_temp;
  int Q_col,Count,size_col,size_col_former;
  double S_m_r_t,S_m_r_t1;
  double C_m,C_1,C_1_m,DELT_S,delt_s;
  int m,L,L_size_col_former_sum;
  double max_S_m_r_t1,min_S_m_r_t1;
  for(t=1;t<=(int)(*t_max);t++)//进入t控制范围
    {
      //把序列分成t个不相交的子序列
      Q=(double *)calloc((int)(*L_Data),sizeof(double));//Q用来存放t个不相交的子序列
      Count=0;
      for(i=0;i<t;i++)
        for(j=0;j<(((int)(*L_Data)-i-1)/t+1);j++)
          Q[Count++]=Data[i+j*t];//把序列分成t个不相交的子序列存放到Q里
      S_m_r_t=0;//记录16个S_m_r_t1的平均值
      DELT_S=0;//记录m的delt_s的平均值
      for(m=2;m<=5;m++) //进入m控制范围
        {
          max_S_m_r_t1=-100;
          min_S_m_r_t1=100;
          for(j=0;j<=3;j++) //进入j控制范围
            {
              S_m_r_t1=0;//记录t个子序列的S_m_r_t1
              L_size_col_former_sum=0;//记录当前子序列以前的所有子序列长度总和
              for(L=0;L<t;L++) //对每个子序列做关联积分计算。
               {
                  
                  if(L==0)
                     size_col_former=0; //用来记录前一个子序列的长度
                  else
                     size_col_former=((int)(*L_Data)-L)/t+1;  //前一个子序列的长度
                     
                  L_size_col_former_sum=L_size_col_former_sum+size_col_former;//记录当前子序列以前的所有子序列长度总和  
                  size_col=((int)(*L_Data)-L-1)/t+1;       //当前的子序列长度
                  
                  Q_temp=(double *)calloc(size_col,sizeof(double));
                  for(j1=0;j1<size_col;j1++)
                    Q_temp[j1]=Q[j1+L_size_col_former_sum]; //当前子序列挪到Q_temp  
                    
                  C_m=correlation_integ(m,1,R[j],size_col,Q_temp);
                  C_1=correlation_integ(1,1,R[j],size_col,Q_temp);
                  
                  C_1_m=1.0;
                  for(j1=0;j1<m;j1++)
                    C_1_m=C_1_m*C_1;
                  S_m_r_t1=S_m_r_t1+(C_m-C_1_m);  
                  free(Q_temp);
                }
                 S_m_r_t1=S_m_r_t1/t;
                 S_m_r_t=S_m_r_t+S_m_r_t1;
                 if(S_m_r_t1>max_S_m_r_t1)
                    max_S_m_r_t1=S_m_r_t1;
                 if(min_S_m_r_t1>S_m_r_t1)
                    min_S_m_r_t1=S_m_r_t1;
            } //j控制范围结束
            delt_s=max_S_m_r_t1-min_S_m_r_t1;
            DELT_S=DELT_S+delt_s;
        } //m控制范围结束
        delt_ave_s[t-1]=DELT_S/4;
        ave_s[t-1]=S_m_r_t/16;
        if(ave_s[t-1]<0)
          S_cor[t-1]=delt_ave_s[t-1]-ave_s[t-1];
        else
          S_cor[t-1]=delt_ave_s[t-1]+ave_s[t-1];
        free(Q);
    }//t控制范围结束
    
 }
 void mexFunction(
     int nlhs,
     mxArray *plhs[],
     int nrhs,
     const mxArray *prhs[]
)
{
	double *x1,*x2,*x3,*x4,*y1,*y2,*y3;
    
	int y1_mrows,y1_ncols;
	int y2_mrows,y2_ncols;
	int y3_mrows,y3_ncols;
	//检查参量的个数
	if (nrhs!=4)
			mexErrMsgTxt("4 inputs required");
  
    x1=mxGetPr(prhs[0]);
	x2=mxGetPr(prhs[1]);
	x3=mxGetPr(prhs[2]);
	x4=mxGetPr(prhs[3]);
	
	  
  y1_mrows=1;
  y1_ncols=(int)(*x3);
  
  y2_mrows=1;
  y2_ncols=y1_ncols;
  
  y3_mrows=1;
  y3_ncols=y1_ncols;
      
    plhs[0]=mxCreateDoubleMatrix(y1_mrows,y1_ncols,mxREAL);
	plhs[1]=mxCreateDoubleMatrix(y2_mrows,y2_ncols,mxREAL);
	plhs[2]=mxCreateDoubleMatrix(y3_mrows,y3_ncols,mxREAL);
      
	y1=mxGetPr(plhs[0]);
	y2=mxGetPr(plhs[1]);
	y3=mxGetPr(plhs[2]);
    DAS(x1,x2,x3,x4,y1,y2,y3);	   			
}