/************************************************************************/
/*                                                                      */
/* Analyse un signal par EMD (Empirical Mode Decomposition)             */
/* i.e. extraction d'IMFs (Intrinsic Mode Function) d'un signal         */
/*                                                                      */
/* Le programme est a compiler a partir de l'invite de commandes Matlab */
/*   >> mex emd.c                                                       */
/*                                                                      */
/* Appel a emd :                                                        */
/*       imfs=emd(x,y,nb_imfs,epsilon);                                 */
/*   ou  imfs=emd(x,y,epsilon,nb_imfs);                                 */
/*                                                                      */
/* Entrees :                                                            */
/*    x : vecteur des abscisses des echantillons du signal a analyser   */
/*    y : vecteur des ordonnees des echantillons du signal a analyser   */
/*    nb_imfs (optionnel) : nombre d'IMFs a extraire du signal          */
/*      -> par defaut, le programme s'arrete lorsqu'il n'y a plus       */
/*         d'IMFs a extraire du signal                                  */
/*    epsilon (optionnel) : parametre de controle de la boucle interne  */
/*                          du programme                                */
/*      -> par défaut, epsilon=0.001 (la moyenne d'un mode est          */
/*         considérée comme nulle si son amplitude crete a crete est    */
/*         inferieure a : epsilon*(amplitude crete a crete du signal y  */
/*         en entree); epsilon appartient a ]0,0.1]                     */
/*                                                                      */
/* Sorties :                                                            */
/*    imfs : matrice contenant les IMFs extraites (une IMF par ligne)   */
/*           ainsi que le residu                                        */
/*                                                                      */
/* Derniere revision : le 11.03.2004                                    */
/*                                                                      */
/************************************************************************/

#include <stdlib.h>
#include "mex.h"
#include "entrees_sorties.h"
#include "min_max.h"
#include "interpolation.h"

#define EPS_DEFAUT 1e-3

double amplitude(double *,int);
int arret(double *,int);

#include "entrees_sorties.c"
#include "min_max.c"
#include "interpolation.c"


/************************************************************************/
/*                                                                      */
/* FONCTION PRINCIPALE                                                  */
/*                                                                      */
/************************************************************************/

void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[]) {

  /* declarations */
  int i,n,choix,nb_imfs,compteur;
  min_max_t mm;
  entree_t entree;
  double *x,*y,*z,*e_max,*e_min,*m,a,eps;
  liste_t liste;

  /* lecture des entrees */
  entree=lire_entrees(nlhs,nrhs,prhs);
  n=entree.n;
  choix=entree.choix;
  nb_imfs=entree.nb_imfs;
  eps=entree.eps;
  x=entree.x;
  y=entree.y;

  /* initialisations */
  mm=initialiser_min_max(n);
  liste=initialiser_liste(n);
  z=(double *)malloc(n*sizeof(double));
  m=(double *)malloc(n*sizeof(double));
  e_max=(double *)malloc(n*sizeof(double));
  e_min=(double *)malloc(n*sizeof(double));
  a=amplitude(y,n);

  /* BOUCLE EXTERNE */

  while ((!choix || (choix && nb_imfs>0)) && !arret(y,n)) {

    /* initialisation */
    for (i=0;i<n;i++) z[i]=y[i];
    for (i=0;i<n;i++) m[i]=y[i];
    compteur=0;

    /* BOUCLE INTERNE */

    while (amplitude(m,n)>=eps*a) {

      /* detection des minima et des maxima */
      min_max(x,z,n,&mm);

      /* interpolation - enveloppe superieure */
      interpolation(e_max,mm.x_max,mm.y_max,mm.n_max,x);

      /* interpolation - enveloppe inferieure */
      interpolation(e_min,mm.x_min,mm.y_min,mm.n_min,x);

      /* calcul de la moyenne */
      for (i=0;i<n;i++) m[i]=(e_max[i]+e_min[i])/2;

      /* mise a jour de z */
      for (i=0;i<n;i++) z[i]=z[i]-m[i];

      /* on compte combien de fois on est passe dans la boucle interne */
      compteur++;

    }

    /* sauvegarde du mode dans la liste si on est     */
    /* passe au moins une fois dans la boucle interne */
    if (compteur) ajouter(&liste,z);

    /* mise a jour de y et nb_imfs */
    for (i=0;i<n;i++) y[i]=y[i]-z[i];
    if (choix) nb_imfs--;

  }

  /* sauvegarde du residu dans la liste */
  ajouter(&liste,y);

  /* ecriture de la sortie dans une matrice */
  ecrire_sorties(liste,plhs);

  /* liberation de la memoire */
  free(y);
  free(m);
  free(e_max);
  free(e_min);
  free(z);
  liberer_liste(liste);
  liberer_min_max(mm);

}


/************************************************************************/
/*                                                                      */
/* AMPLITUDE D'UN VECTEUR                                               */
/*                                                                      */
/* Fonction appelee par la fonction principale pour evaluer le critere  */
/* d'arret de la boucle interne.                                        */
/* Elle calcule la difference entre le maximum global et le minimum     */
/* global d'un vecteur.                                                 */
/*                                                                      */
/************************************************************************/

double amplitude(double *u,int n) {
  double max,min;
  int i;
  min=max=u[0];
  for (i=1;i<n;i++) {
    if (u[i]>max) max=u[i];
    if (u[i]<min) min=u[i];
  }
  return max-min;
}


/************************************************************************/
/*                                                                      */
/* TEST D'ARRET DE LA BOUCLE EXTERNE                                    */
/*                                                                      */
/* Fonction appelee par la fonction principale pour evaluer le critere  */
/* d'arret de la boucle externe.                                        */
/* Elle teste s'il reste au moins une periode dans le signal.           */
/*                                                                      */
/************************************************************************/

int arret(double y[],int n) {
  int cour,min=0,max=0;
  cour=1;

  /* comptage des minima et des maxima */
  while (cour<(n-1)) {
    if (y[cour]<=y[cour-1] && y[cour]<=y[cour+1]) min++;
    if (y[cour]>=y[cour-1] && y[cour]>=y[cour+1]) max++;
    cour++;
    while (y[cour]==y[cour-1] && cour<(n-1)) cour++;
  }
  if (n>=3 && y[0]==y[1]) {
    if (y[2]>=y[1]) min--; /* on a compte, a tort, y[1] comme un minimum */
    if (y[2]<=y[1]) max--; /* on a compte, a tort, y[1] comme un maximum */
  }
  if (n>=3 && y[n-1]==y[n-2]) {
    if (y[n-3]>=y[n-2]) min--; /* on a compte, a tort, y[n-2] (*) */
			       /* comme un minimum                */
    if (y[n-3]<=y[n-2]) max--; /* on a compte, a tort, y[n-2] (*) */
                               /* comme un maximum                */
  }

  /* (*) ou un autre avant lui, ayant la meme valeur */

  /* s'il y a moins de trois extrema, on s'arrete */
  if ((min>=2 && max>=1) || (max>=2 && min>=1)) return 0;
  else return 1;
}
