#include "mex.h"
#include "matlab.h"
#include "math.h"
#define LOG(a) ((a)>1e-20 ? log((a)):log(1e-20))
#include <afxwin.h> // MFC core and standard components
#include <afxext.h> // MFC extensions
#include <afxdisp.h> // MFC Automation classes
#include <afxdtctl.h> // MFC support for Internet Explorer 4 Common Controls
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
if((nrhs!=3)||(nlhs!=1))
{
mexPrintf("********************--------correl_dim----------*******************\n");
mexPrintf("本函数计算相关积供后续程序计算关联维\n");
mexPrintf("c=correl_dim(data,n_eps,stepsize)\n");
mexPrintf("输入：data-------待处理数据 相点坐标维数*点数；\n");
mexPrintf(" n_eps---搜寻径点数 1*1; \n");
mexPrintf(" stepsize---搜寻径步长 1*1; eps[l]=maxdist*2.0/pow(buchang,n_eps-l);\n");
mexPrintf("输出：c----------相关积结果 搜寻径点数*1;\n");
mexPrintf("********************--------correl_dim----------*******************\n");
mexErrMsgTxt("你想怎么用亚");
}
int Mr0,Nr0;
Mr0=mxGetM(prhs[0]);
Nr0=mxGetN(prhs[0]);
if(Nr0<102)
mexErrMsgTxt("数据点数过少，不宜计算。");//延迟100个数据点，对于特征数据这个是不必要的。
int n_eps=mxGetScalar(prhs[1]);
if(n_eps<1)
mexErrMsgTxt("n_eps<1");
double buchang=mxGetScalar(prhs[2]);
double *y,*jieguo;
plhs[0]=mxCreateDoubleMatrix(n_eps,1,mxREAL);
y =mxGetPr(prhs[0]);
jieguo=mxGetPr(plhs[0]);
double *eps=(double*)calloc(n_eps,sizeof(double));
for(int l=0;l<n_eps;l++)
{
eps[l]=0;
jieguo[l]=0;
}
double maxdist=0;
for(int j=0;j<Nr0;j++)
{
double dist=0;
double *a=y;
double *b=y+j*Mr0;
for(int mm=0;mm<Mr0;mm++)
dist=dist+(a[mm]-b[mm])*(a[mm]-b[mm]);
if(maxdist<dist)
maxdist=dist;
}
if(maxdist<=0)
mexErrMsgTxt("数都一样");
for(l=0;l<n_eps;l++)
eps[l]=maxdist*2.0/pow(buchang,n_eps-l);
for(j=0;j<Nr0-100;j++)
{
for(int k=j+100;k<Nr0;k++)
{
double dist=0;
double *a=y+j*Mr0;
double *b=y+k*Mr0;
for(int mm=0;mm<Mr0;mm++)
dist=dist+(a[mm]-b[mm])*(a[mm]-b[mm]);
for(int l=0;l<n_eps;l++)
{
if(dist<=eps[l])
jieguo[l]+=1;
}
}
}
for(l=0;l<n_eps;l++)
{
jieguo[l]=LOG(jieguo[l])/(0.5*LOG(buchang));
}
free(eps);
}
