g-p.txt

G-P算法程序

该程序转载于《分形与混沌——blog》
再次声明是转载!!!!
为了自己的学习，记录别人的好的东东，！
 
#include "mex.h"
#include "math.h"
/*****************************************************************
*  data:the time series
*  r: the r series
*  N: the length of the time series
*  delay : the delay of the time series
*  rlen: the length of r series
*  ln_C: the correlation dimention series
******************************************************************/
void G_P(double data[], double r[], int N, int m,
   int rlen, int delay, double ln_C[],double ln_r[])
{
 double distance;
 double max=0;
 int cr=0;
 int i,j,k,l;
 double Cr=0;
 int M=N-(m-1)*delay;
 for (l=0; l<rlen; l++)
 {
     for (i=0; i<M-1; i++)
  {
      for (j=i+1; j<M; j++)
   {
       for (k=0; k<m; k++)
    {
    distance = fabs(data[i+k*delay]-data[j+k*delay]);
    if (distance>max)
     max = distance;
    }
       if (r[l]>=max)
        cr++;
       max = 0;
   }
  }
  Cr = cr*2.0/(M*(M-1));
  cr = 0;
        ln_C[l] = log(Cr);
        ln_r[l] = log(r[l]);  
 }
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
 double *data;
 double *r;
 int m;
 int delay;
 double *ln_C;
 double *ln_r;
 int datalen,rlen;
 if (nrhs!=4)
  mexErrMsgTxt("Four inputs required!"); 
 else if (nlhs>2)
  mexErrMsgTxt("Too many output arguments!");
 datalen = mxGetN(prhs[0]);
 rlen = mxGetN(prhs[1]);
 m = mxGetScalar(prhs[2]);
 delay = mxGetScalar(prhs[3]);
    //define the two output arguments
 plhs[0] = mxCreateDoubleMatrix(1,rlen,mxREAL);
 plhs[1] = mxCreateDoubleMatrix(1,rlen,mxREAL);
 data = mxGetPr(prhs[0]);
 r = mxGetPr(prhs[1]);
 ln_C = mxGetPr(plhs[0]);
 ln_r = mxGetPr(plhs[1]);
 G_P(data,r,datalen,m,rlen,delay,ln_C,ln_r);
}