该程序转载于《分形与混沌——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);
}
评论0