lyap_exp_wolf.cpp

本函数计算距离增长结果供拟合lyapunov指数 详情见说明

#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("********************--------lyap_exp----------*******************\n");
		mexPrintf("本函数计算距离增长结果供拟合lyapunov指数\n");
		mexPrintf("c=lyap_exp_wolf(data,yanchi,t)\n");
		mexPrintf("输入：0---data-------待处理数据 相点坐标维数*点数；\n");
		mexPrintf("      1---yanchi-----延迟点数   1*1;\n");
		mexPrintf("      2---t----------采样时间   1*1;\n");
		mexPrintf("输出：0---c----------结果       yanchi*n_eps;\n");
		mexPrintf("基于wolf轨线追踪方法\n");
		mexPrintf("********************--------lyap_exp----------*******************\n");
		mexErrMsgTxt("你想怎么用亚");
	}
	    
	int Mr0,Nr0;
	Mr0=mxGetM(prhs[0]);
	Nr0=mxGetN(prhs[0]);
	int yanchi=mxGetScalar(prhs[1]);
	if(Nr0<100+yanchi)
		mexErrMsgTxt("数据点数过少，不宜计算。");//延迟100个数据点，对于特征数据这个是不必要的。

	double t=mxGetScalar(prhs[2]);
	
	double *y,*jieguo;
	plhs[0]=mxCreateDoubleMatrix(yanchi,1,mxREAL);
	y     =mxGetPr(prhs[0]);
	jieguo=mxGetPr(plhs[0]);	
	
	for(int j=0;j<yanchi;j++)
		jieguo[j]=0;	

	int yongdaodianshu=0;
	for(j=0;j<Nr0-yanchi;j++)
	{
		int mink=0;
		double mindist=0;
		for(int mm=0;mm<Mr0;mm++)
			mindist=mindist+(y[j*Mr0+mm]-y[mm])*(y[j*Mr0+mm]-y[mm]);
		
		for(int k=0;k<Nr0-yanchi;k++)
		{
			if (abs(k-j)<10) 
				continue;
			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]);
			if(dist<mindist)
			{
				mindist=dist;   
				mink=k;
			}
		}
		
		while(mindist<1e-10)
		{
			mink=k++;
			if (mink==Nr0-yanchi) 
				break;
			double *a=y+j*Mr0;
			double *b=y+mink*Mr0;
			for(int mm=0;mm<Mr0;mm++)
				mindist=mindist+(a[mm]-b[mm])*(a[mm]-b[mm]);
		}
		if(mindist>=1e-10)
		{
			for(int yan=0;yan<yanchi;yan++)
			{
				double *a_yan=y+(j+yan)*Mr0;
				double *b_yan=y+(mink+yan)*Mr0;
				double dist=0;
				for(int mm=0;mm<Mr0;mm++)
					dist=dist+(a_yan[mm]-b_yan[mm])*(a_yan[mm]-b_yan[mm]);
				jieguo[yan]+=LOG(dist/mindist);
			}
			yongdaodianshu++;
		}
	}
	for(int qq=0;qq<yanchi;qq++)
		jieguo[qq]/=(t*2.0*yongdaodianshu);
}