multi_correl_dim.cpp

本函数计算相关积供后续程序计算计算多重关联维 详情见程序内说明

#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!=4)||(nlhs!=1))
	{
		mexPrintf("********************--------multi_correl_dim----------*******************\n");
		mexPrintf("本函数计算相关积供后续程序计算多重关联维\n");
		mexPrintf("c=multi_correl_dim(data,n_eps,stepsize,q)\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("      q----------多重指标;  1*qm;");
		mexPrintf("输出：c----------相关积结果 搜寻径点数*qm;\n");
		mexPrintf("********************--------multi_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]);
	int Mr3=mxGetN(prhs[3]);
	double *q=mxGetPr(prhs[3]);
		
	double *y,*jieguo;
	plhs[0]=mxCreateDoubleMatrix(n_eps,Mr3,mxREAL);
	y     =mxGetPr(prhs[0]);
	jieguo=mxGetPr(plhs[0]);	
	double *eps=(double*)calloc(n_eps,sizeof(double));	
	double *zhongzhuan=(double*)calloc(n_eps,sizeof(double));
	
	for(int l=0;l<n_eps;l++)
		eps[l]=0;
	for(l=0;l<n_eps*Mr3;l++)
		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=11;j<Nr0-11;j++)
	{
		for(int k=0;k<j-10;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])
					zhongzhuan[l]+=1;	
			}
		}
		
		for(k=j+10;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])
					zhongzhuan[l]+=1;	
			}
		}

		for(int qq=0;qq<Mr3;qq++)
		{
			for(int jj=0;jj<n_eps;jj++)
			{
				if(q[qq]==1)
					q[qq]=2;
				jieguo[qq*n_eps+jj]+=pow(zhongzhuan[jj],q[qq]-1);
			}
		}
	}
	for(int qq=0;qq<Mr3;qq++)
	{
		for(int jj=0;jj<n_eps;jj++)
		{
			jieguo[qq*n_eps+jj]=2.0*LOG(jieguo[qq*n_eps+jj])/((q[qq]-1)*LOG(buchang));
		}
	}
	free(eps);
	free(zhongzhuan);
}