cio.c

经验模式分解源码 学习希尔波特黄变换的参考程序 源于Internet网络

/*
* G. Rilling, last modification: 3.2007
* [e-mail]
*
* code based on a student project by T. Boustane and G. Quellec, 11.03.2004
* supervised by P. Chainais (ISIMA - LIMOS - Universite Blaise Pascal - Clermont II
* email : [e-mail]).
*/

/************************************************************************/
/*                                                                      */
/* GET INPUT DATA                                                       */
/*                                                                      */
/************************************************************************/

input_t get_input(int nlhs,int nrhs,const mxArray *prhs[]) {
  input_t input;
  int n,i;
  double *x,*ry_temp,*iy_temp,*third,fourth,fifth;
  COMPLEX_T *y;
  
  input.stop_params.threshold = DEFAULT_THRESHOLD;
  input.stop_params.tolerance = DEFAULT_TOLERANCE;
  input.allocated_x=0;
  #ifdef _ALT_MEXERRMSGTXT_
  input.error_flag=0;
  #endif
  input.max_imfs=0;
  input.nbphases=DEFAULT_NBPHASES;
  
  /* argument checking*/
  if (nrhs>5)
    mexErrMsgTxt("Too many arguments");
  if (nrhs<2)
    mexErrMsgTxt("Not enough arguments");
  if (nlhs>2)
    mexErrMsgTxt("Too many output arguments");
  if (!mxIsEmpty(prhs[0]))
    if (!mxIsNumeric(prhs[0]) || mxIsComplex(prhs[0]) ||
    mxIsSparse(prhs[0]) || !mxIsDouble(prhs[0]) ||
    (mxGetNumberOfDimensions(prhs[0]) > 2))
      mexErrMsgTxt("X must be either empty or a double precision real vector.");
  
  if (!mxIsNumeric(prhs[1]) || !mxIsComplex(prhs[1]) ||
  mxIsSparse(prhs[1]) || !mxIsDouble(prhs[1]) ||/* length of vector x */
  (mxGetNumberOfDimensions(prhs[1]) > 2))
    mexErrMsgTxt("Y must be a double precision complex vector.");
  
  /* input reading: x and y */
  n=GREATER(mxGetN(prhs[1]),mxGetM(prhs[1])); /* length of vector x */
  if (mxIsEmpty(prhs[0])) {
    input.allocated_x = 1;
    x = (double *)malloc(n*sizeof(double));
    for(i=0;i<n;i++) x[i] = i;
  }
  else
    x=mxGetPr(prhs[0]);
  ry_temp=mxGetPr(prhs[1]);
  iy_temp=mxGetPi(prhs[1]);
  
  /* third argument */
  if (nrhs>=3) {
    if(!mxIsEmpty(prhs[2])) {
      if (!mxIsNumeric(prhs[2]) || mxIsComplex(prhs[2]) || mxIsSparse(prhs[2])
      || !mxIsDouble(prhs[2]) || (mxGetN(prhs[2])!=1 && mxGetM(prhs[2])!=1))
        mexPrintf("STOP must be a real vector of 1 or 2 elements");
      i = GREATER(mxGetN(prhs[2]),mxGetM(prhs[2]));
      if (i>2)
        mexErrMsgTxt("STOP must be a vector of 1 or 2 elements");
      third=mxGetPr(prhs[2]);
      switch (i) {
        case 1 : {
          input.stop_params.threshold=*third;
          break;
        }
        case 2 : {
          input.stop_params.threshold=third[0];
          input.stop_params.tolerance=third[1];
        }
      }
      /* input checking */
      if (input.stop_params.threshold <= 0)
        mexErrMsgTxt("threshold must be a positive number");
      if (input.stop_params.threshold >= 1)
        mexWarnMsgTxt("threshold should be lower than 1");
      if (input.stop_params.tolerance < 0 || input.stop_params.tolerance >= 1)
        mexErrMsgTxt("tolerance must be a real number in [O,1]");
    }
  }
  
  
  /* fourth argument */
  if (nrhs>=4) {
    if (!mxIsEmpty(prhs[3])) { /* if empty -> do nothing */
      if (!mxIsNumeric(prhs[3]) || mxIsComplex(prhs[3]) || mxIsSparse(prhs[3])
      || !mxIsDouble(prhs[3]) || mxGetN(prhs[3])!=1 || mxGetM(prhs[3])!=1)
        mexErrMsgTxt("NB_IMFS must be a positive integer");
      fourth=*mxGetPr(prhs[3]);
      if ((unsigned int)fourth != fourth)
        mexErrMsgTxt("NB_IMFS must be a positive integer");
      input.max_imfs=(int)fourth;
    }
  }
  
  /* fifth argument */
  if (nrhs==5) {
    if(!mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) || mxIsSparse(prhs[4])
    || mxGetN(prhs[4])!=1 || mxGetM(prhs[4])!=1)
      mexErrMsgTxt("NBPHASES must be a positive integer");
    fifth=*mxGetPr(prhs[4]);
    if ((int)fifth != fifth)
      mexErrMsgTxt("NBPHASES must be a positive integer");
    input.nbphases = (int)fifth;
  }
  
  
  /* more input checking */
  if (!input.allocated_x && SMALLER(mxGetN(prhs[0]),mxGetM(prhs[0]))!=1 ||
  SMALLER(mxGetN(prhs[1]),mxGetM(prhs[1]))!=1)
    mexErrMsgTxt("X and Y must be vectors");
  if (GREATER(mxGetN(prhs[1]),mxGetM(prhs[1]))!=n)
    mexErrMsgTxt("X and Y must have the same length");
  i=1;
  while (i<n && x[i]>x[i-1]) i++;
  if (i<n) mexErrMsgTxt("Values in X must be non decreasing");
  
  /* copy vector y to avoid erasing input data */
  y=(COMPLEX_T *)malloc(n*sizeof(COMPLEX_T));
  #ifdef C99_OK
  for (i=0;i<n;i++) y[i]=ry_temp[i]+I*iy_temp[i];
  #else
  for (i=0;i<n;i++) {
    y[i].r=ry_temp[i];
    y[i].i=iy_temp[i];
  }
  #endif
  
  input.n=n;
  input.x=x;
  input.y=y;
  return input;
}


/************************************************************************/
/*                                                                      */
/* INITIALIZATION OF THE LIST                                           */
/*                                                                      */
/************************************************************************/

imf_list_t init_imf_list(int n) {
  imf_list_t list;
  list.first=NULL;
  list.last=NULL;
  list.n=n;
  list.m=0;
  return list;
}


/************************************************************************/
/*                                                                      */
/* ADD AN IMF TO THE LIST                                               */
/*                                                                      */
/************************************************************************/

void add_imf(imf_list_t *list,COMPLEX_T *p,int nb_it) {
  COMPLEX_T *v=(COMPLEX_T *)malloc(list->n*sizeof(COMPLEX_T));
  int i;
  imf_t *mode=(imf_t *)malloc(sizeof(imf_t));
  for (i=0;i<list->n;i++) v[i]=p[i];
  mode->pointer=v;
  mode->nb_iterations=nb_it;
  mode->next=NULL;
  if (!list->first) {
    list->first=mode;
  } else {
    (list->last)->next=mode;
  }
  list->last=mode;
  list->m++;
}


/************************************************************************/
/*                                                                      */
/* FREE MEMORY ALLOCATED FOR THE LIST                                   */
/*                                                                      */
/************************************************************************/

void free_imf_list(imf_list_t list) {
  imf_t *current=list.first, *previous;
  while (current) {
    previous=current;
    current=current->next;
    free(previous->pointer);
    free(previous);
  }
}


/************************************************************************/
/*                                                                      */
/* OUTPUT INTO MATLAB ARRAY                                             */
/*                                                                      */
/************************************************************************/

void write_output(imf_list_t list,mxArray *plhs[]) {
  double *rout,*iout,*out2;
  imf_t *current;
  int i=0,j,m=list.m,n=list.n;
  plhs[0]=mxCreateDoubleMatrix(m,n,mxCOMPLEX);
  rout=mxGetPr(plhs[0]);
  iout=mxGetPi(plhs[0]);
  plhs[1]=mxCreateDoubleMatrix(1,m-1,mxCOMPLEX);
  out2=mxGetPr(plhs[1]);
  for (current=list.first;current;current=current->next) {
    for (j=0;j<n;j++) {
      *(rout+j*m+i)=CREAL(current->pointer[j]);
      *(iout+j*m+i)=CIMAG(current->pointer[j]);
    }
    if (i<m-1) *(out2+i)=current->nb_iterations;
    i++;
  }
}