lyapunov_rosenstein_c.c

Hopfield神经网络算法源程序：这程序是《神经网络模式识别及其实现》第九章<神经联想记忆和Hopfield网络>附的程序

#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]        // 范数类型 norm = 0,1,2时，分别对应无穷范数、1范数和2范数
#define P prhs[5]           // 限制短暂分离，大于序列平均周期(不考虑该因素时 p = 1) 
#define S prhs[6]           // 采样频率


// 定义输出参数
#define LYAP plhs[0]

// 声明 C 运算函数 (该函数名不能和本文件名重名)
void LYAPUNOV_EXPONENT();
double NORM_VECTOR();
int MIN();

void mexFunction (int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])	
{
    double *px,*pc;
    int m,t,s,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);     // 范数类型 norm = 0,1,2时，分别对应无穷范数、1范数和2范数
    p = (int) mxGetScalar(P);           // 限制短暂分离，大于序列平均周期(不考虑该因素时 p = 1)
    s = (int) mxGetScalar(S);           // 采样频率

    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;
    }
    // 调用 C 运算函数 (该函数名不能和本文件名重名)
    LYAPUNOV_EXPONENT(px,N,m,length,norm,p,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))*s;
    	}

    }
}


//--------------------------------------------------------------------------
//           计算lyapunov指数曲线。特别注意：这里输入参数为重构矩阵
//--------------------------------------------------------------------------
void LYAPUNOV_EXPONENT(double *pxn,     // 重构矩阵 
                          int xn_cols,  // 相空间点数
                          int m,        // 嵌入维数 
                          int length,   // 最长演化时间
                          int norm,     // 范数类型 norm = 0,1,2时，分别对应无穷范数、1范数和2范数 
                          int p,        // 限制短暂分离，大于序列平均周期(不考虑该因素时 p = 1)
                          int *evolution_sum,
                          double *pc)        
{
    double d_ij,*pd_vector,temp,c,d,min_distance,log_d_ij;
    int i,j1,j2,j3,distance_index,max_evolution,evolution_fact;
   
    pd_vector = malloc(m*sizeof(double));               // 重构向量 m*1
    temp = 0;
    distance_index=0;
    
    for (j1=0; j1<xn_cols; j1++)
    {   
    	min_distance=100;
    	for (j2=0; j2<xn_cols; j2++)
        {   
        	if ((abs(j1-j2)>p)&&(abs(j1-j2)<30*p))
        	{              	
        		for (i=0;i<m;i++)
        		{
        			d = *(pxn+i*xn_cols+j1) - *(pxn+i*xn_cols+j2);    // 第j1列与j2列对应元素相减
        			*(pd_vector+i) = d;
        		}
        		// 计算数组的范数　norm = 0,1,2时，分别对应无穷范数、1范数和2范数
                        d_ij = NORM_VECTOR(pd_vector,m,norm);
                        if (d_ij<min_distance)
                        {
                        	min_distance=d_ij;
                        	distance_index=j2;
                        }
                }
        }
        max_evolution=MIN((xn_cols-j1),(xn_cols-distance_index));           //计算点j的最大演化时间步长i，这里的演化下一步中进行了改正
	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+distance_index+j3);    // 第j1列与j2列对应元素相减
        		*(pd_vector+i) = d;
        	}
                // 计算数组的范数　norm = 0,1,2时，分别对应无穷范数、1范数和2范数
                d_ij = NORM_VECTOR(pd_vector,m,norm);
                if (d_ij<0.000001)
                                d_ij=0.000001;
                log_d_ij=log(d_ij);
                *(pc+j3)=*(pc+j3)+log_d_ij;
                *(evolution_sum+j3)=*(evolution_sum+j3)+1;
        }
    }
}


//---------------------------------------------------------------------------
// 计算数组的范数　norm = 0,1,2时，分别对应无穷范数、1范数和2范数
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:     // 1 - 范数
            value = 0;
            for (i=1; i<len; i++)
            {   
                value = value + abs(*(p_vector+i));            
            }
            break;
            
        case 2:     // 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;
}