/** G. Rilling, last modification: 3.2007* gabriel.rilling@ens-lyon.fr** 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 : pchainai@isima.fr).*//************************************************************************//*                                                                      *//* 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++;  }}