cao1.c

cao_m.m----主程序, cao法求最小嵌入维程序 reconstitution.m---------子程序,相重构程序

#include <math.h>
#include "mex.h"
#include "stdio.h"
#include "stdlib.h"
#include "matrix.h"

//---------------------------------------------------------------------------
// 计算绝对值
double ABS(double x)
{
    double y;
    if (x>=0)
        y = x;
    else
        y = -x;
    return y;
}
//---------------------------------------------------------------------------

// 计算数组最大值
double MAX_VECTOR(double *p_vector,
                  int len_vector)
{
    int i;    
    double max_value = *p_vector;
    for (i=0; i<len_vector; i++)
    {   if (*(p_vector+i)>max_value)
	{
		max_value = *(p_vector+i);
	}
    }
    return max_value;
}




//求最短距离及向量
void min_dist(double *pdist,int n,int *place,double *pmindist,int i)
{  int k,j;    
   double min_value ;
   
   for (k=0;k<n;k++)
   {
	   if (*(pdist+k)!=0) 
	   {min_value=*(pdist+k);j=k+1;
	   break;}
   }

  for (k=0; k<n; k++)
  {   
	  if ((*(pdist+k)!=0)&&(*(pdist+k)<min_value))
	  {	 min_value = *(pdist+k); j=k+1; }
  }
  
  *(place+i)=j;
  *(pmindist+i)=min_value;

}



//-----------------------------------------------------------------------------
// 计算向量间的距离（最大范数）
double distance(double *Y1,double *Y2,int m1,int tau)
{  int i;
  double d,*pd,d_ij;
  pd = (double*)malloc(m1*sizeof(double));       // 声明一个长度为 m 的 double 型数组
	for (i=1;i<=m1;i++)
	{
			d=*(Y1+(i-1)*tau)-*(Y2+(i-1)*tau);
			*(pd+i-1)=ABS(d);
	}
	d_ij=MAX_VECTOR(pd,m1);
    free(pd);
	return d_ij;
}


//---------------------------------------------------------------------------

// 定义 C 运算函数
void Cao_FUNCTION( double *pdata,      // 时间序列（列向量）
					int m_min,          // 最小嵌入维数       
					int m_max,          // 最大嵌入维数 
					int tau,            // 时间延迟
					int N ,  		     // 序列长度
					double *pE1,        //用于存放E1
					double *pE2)        //用于存放E2
					
{   int *place,i,j,m;
    double *pdist,*pmindist,*pE,*pEn;

   	pE = (double*)malloc((m_max-m_min+1)*sizeof(double)); 
	pEn = (double*)malloc((m_max-m_min+1)*sizeof(double)); 

 for (m=m_min;m<=m_max;m++)
 {    
     double ad=0,neard, dd=0;
      int n=N-m*tau,d;
      place = (int*) malloc(n*sizeof(int));    //用于存放重构后向量的标号
	  pmindist = (double*)malloc(n*sizeof(double)); 
	 // printf("共有 %d 步,正在计算第 %d 步!\n",m_max,m);
	 	  
    for (i=0;i<n;i++)         //计算第i个向量的距离最近的向量
	{  
		pdist = (double*)malloc(n*sizeof(double));   //存放第i个与每个向量的距离

       for (j=0;j<n;j++)
	   { *(pdist+j)=distance(pdata+i,pdata+j,m,tau);	}
	  
	    min_dist(pdist,n,place,pmindist,i);   //求第i个向量的最近距离及对应向量标号

      //--------------计算a(i,m)-------------------
      d=*(place+i);
      neard=*(pmindist+i);
      ad=ad+distance(pdata+i,pdata+d-1,m+1,tau)/neard;
	  dd=dd+ABS(*(pdata+i+m*tau)-*(pdata+d-1+m*tau));//求E*(d)
      free(pdist);
	}
	
      *(pE+m-m_min)=ad/n;   //求E(d)
      *(pEn+m-m_min)=dd/n; 	 //求E*(d)
      
      free(place);free(pmindist);
 }
//----------求E1(d)&E2(d)-------------------
  for (i=1;i<=m_max-m_min;i++)
  {
	   *(pE1+i-1)=*(pE+i)/(*(pE+i-1));
	   *(pE2+i-1)=*(pEn+i)/(*(pEn+i-1));
	   
  }
  free(pE);free(pEn);
}


//---------------------------------------------------------------------------
void 
mexFunction (int nlhs, mxArray *plhs[],			// 输出参数个数，及输出参数数组
			 int nrhs, const mxArray *prhs[])	// 输入参数个数，及输入参数数组
{
      double *pdata,*pE1,*pE2;
      int m_min,m_max,tau,N,nRows;

	  if (nrhs!=4) mexErrMsgTxt("需要4个参数!");  //检查输入参数的个数

    // 取得输入参数
    pdata = mxGetPr(prhs[0]);      // 时间序列（列向量）      
    m_min = (int) *mxGetPr(prhs[1]);            // 最小嵌入维数 
	m_max = (int) *mxGetPr(prhs[2]);            // 最大嵌入维数 
    tau = (int) *mxGetPr(prhs[3]);      // 时间延迟   
	N = mxGetM(prhs[0]);              // 序列长度
   
	
    // 为输出变量分配内存空间
	nRows=m_max-m_min;
	 plhs[0]= mxCreateDoubleMatrix(nRows,1,mxREAL); //用于存放E1
	 plhs[1]= mxCreateDoubleMatrix(nRows,1,mxREAL); //用于存放E2
	
	// 取得输出参数指针
	pE1 = mxGetPr(plhs[0]);
    pE2 = mxGetPr(plhs[1]);

    // 调用 C 运算函数 (该函数名不能和本文件名重名)
     Cao_FUNCTION(pdata,m_min,m_max,tau,N,pE1,pE2);    
	
}