lyapunov_kantz_c.c

混沌kantz算法的mex实现kantz_lyapunov_c

#include <math.h>
#include "mex.h"
#include "stdio.h"
#include "stdlib.h"
#include "matrix.h"


#define X prhs[0]           
#define T prhs[1]           
#define M prhs[2]           
#define EVLU prhs[3]        
#define NORM prhs[4]        
#define P prhs[5]           
#define FS prhs[6]           
#define R prhs[7]           

#define LYAP plhs[0]

void LYAPUNOV_EXPONENT();
double NORM_VECTOR();
int MIN();

void mexFunction (int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])	
{
    double *px,*pc,radius;
    int m,t,fs,N,p,norm,length,m_index;
    int *evolution_sum;
    
    px = mxGetPr(X);                
    N = mxGetM(X);                      
    t = (int) mxGetScalar(T);            
    m = (int) mxGetScalar(M);                
    length = (int) mxGetScalar(EVLU);     
    norm = (int) mxGetScalar(NORM);   
    p = (int) mxGetScalar(P);         
    fs = (int) mxGetScalar(FS);          
    radius = mxGetScalar(R);        
    
    evolution_sum=(int*)malloc(length*sizeof(int));
    
   LYAP = mxCreateDoubleMatrix(1,length,mxREAL); 
   pc = mxGetPr(LYAP);
	
    for (m_index=0;m_index<length;m_index++)
    {
    	*(evolution_sum+m_index)=0;
    	*(pc+m_index)=0;
    }

    LYAPUNOV_EXPONENT(px,N,m,length,norm,p,radius,evolution_sum,pc);    
    for (m_index=0;m_index<length;m_index++)
    {
	if(*(evolution_sum+m_index)>0)
    	{
    		*(pc+m_index)=(*(pc+m_index))/(*(evolution_sum+m_index))*fs;
    		//mexPrintf("we have counted %d numbers of dots %d.\n", *(evolution_sum+m_index),*(pc+m_index));
    	}
    }
}

void LYAPUNOV_EXPONENT(double *pxn,     
                          int xn_cols,  
                          int m,        
                          int length,   
                          int norm,    
                          int p,    
                          double radius,  
                          int *evolution_sum,
                          double *pc)        
{
    double d,d_ij,*pd_vector,log_d_ij,d_ij_revolution;
    int i,j1,j2,j3,max_evolution,evolution_fact;

    pd_vector = malloc(m*sizeof(double));      
      
    for (j1=0; j1<xn_cols; j1++)
    {   

    	for (j2=0; j2<xn_cols; j2++)
        {   
        	if ((abs(j1-j2)>p)&&(abs(j1-j2)<10*p))
        	{           	
        		for (i=0;i<m;i++)
        		{
        			d = *(pxn+i*xn_cols+j1) - *(pxn+i*xn_cols+j2);   
        			*(pd_vector+i) = d;
        		}

                        d_ij = NORM_VECTOR(pd_vector,m,norm);
                        if (d_ij<radius)
                        {
                        	max_evolution=MIN((xn_cols-j1),(xn_cols-j2));          
	                        evolution_fact=MIN(length,max_evolution); 
	
	                        for (j3=0;j3<evolution_fact;j3++)   
	                        {
	                        	for (i=0;i<m;i++)
	                        	{
	                        		d = *(pxn+i*xn_cols+j1+j3) - *(pxn+i*xn_cols+j2+j3);    
        		                        *(pd_vector+i) = d;
        	                        }

                                        d_ij_revolution=NORM_VECTOR(pd_vector,m,norm);     
                                        if (d_ij_revolution<0.000001)
                                               d_ij_revolution=0.000001;
                                        log_d_ij=log(d_ij_revolution);
                                        *(pc+j3)=*(pc+j3)+log_d_ij;
                                        *(evolution_sum+j3)=*(evolution_sum+j3)+1;      
                               }
                        }
                  }
            }  
       //mexPrintf("we have counted %d numbers of dots.\n", j1);
       }

       free(pd_vector);
}


double NORM_VECTOR(double *p_vector,int len,int norm)
{
    int i;    
    double value,tmp;
    
    switch (norm)
    {
        case 0:     
            value = abs(*(p_vector));
            for (i=2; i<len; i++)
            {   
                tmp = abs(*(p_vector+i));            
                if (tmp>value)
                {
                    value = tmp;
                }
            }
            break;
   
        case 1:     
            value = 0;
            for (i=1; i<len; i++)
            {   
                value = value + abs(*(p_vector+i));            
            }
            break;
            
        case 2:     
            value = 0;
            for (i=1; i<len; i++)
            {   
                value = value + (*(p_vector+i))*(*(p_vector+i));
            }
            value = sqrt(value);
            break;
    }
    return value;
}

int MIN(int x,int y)
{
    
    if (x<=y)
        return x;
    else
        return y;
}