#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <cstring>

#include "header.h"
#include "PlayWav.h"
#include "FFT.h"
#include "Functions.h"

using namespace std;

//Fichier WAV
struct s_wSamples WAV_wSamples;
struct s_wVector MFCC_wVector1;
struct s_wVector MFCC_wVector2;
FILE *WAV_FILE;
float DTW_value, DTW_tmp;

const char tmp_file [8] = "tmp.mfc";

//Programme principal
int main(int argc, char *argv[])
{
  char *input_file = NULL;
  char *output_file = NULL;
  char *save_file = NULL;
  int error = -1;
  int sError = -1;
  bool sPlay = false;
  bool sMFCC = false;
  bool sDTW = false;
  bool sDICO = false;
  float *DataOut = NULL;
  char buffer [256], buffer_file [256];
  
  //Efface écran affiche texte
  //system ("CLS");
  printf ("RobotVoice v0.21\n");
  printf ("This program run with the library FMOD\n");
  printf ("Copyright 2004-2005 by ANNEHEIM Geoffrey, MONCHAMPS Sebastien,\n");
  printf ("                       DE LA BROISE Jeremi and DEPONT Jacques.\n\n");
  printf ("Command:\n");
  printf ("-p = Play PCM file.\n");
  printf ("-s = Make MFCC file.\n");
  printf ("-c = Compare two MFCC files.\n");
  printf ("-d = Compare WAV file with a dictionnary.\n\n");
  printf ("Usage: MFCC file or PCM file or Text file [-p][-s][-c][-d]\n\n");
  
   //Vide structure
  memset (&WAV_wSamples, 0, sizeof (s_wSamples));
  
  //Charge la librairie FMOD
  if (!PlayWav_Initialize (SAMPLE_RATE, MAX_CHANNELS)) {
    printf ("Error during initialization of FMOD.\n");
    system ("PAUSE");
    return 4;}            
  
  //Parcours arguments
  for (int i=1; i<argc; i++)
  {
    //==================
    //PLAY FILE     [-p]
    //==================
    if (!strcmp (argv [i], "-p") && sPlay==false)
    {
      sPlay = true;
      continue;
    }
    //===================
    //MAKE MFCC FILE [-s]
    //===================
    else if (!strcmp (argv [i],"-s") && sMFCC==false)
    {
      sMFCC = true;
      continue; 
    }
    //==================
    //COMPARE (DTW) [-c]
    //==================
    else if (!strcmp (argv [i],"-c") && sDTW==false)
    {
      sDTW = true;
      continue;    
    } 
    //=====================
    //COMPARE WITH DIC [-d]
    //=====================
    else if (!strcmp (argv [i],"-d") && sDICO==false)
    {
      sDICO = true;
      continue;    
    }       
    //Argument incorrect
    else 
    {
      if (!input_file) input_file = argv [i];
      else if (!output_file) output_file = argv [i];
      else {
        printf ("Error: Incorrect argument.\n");
        sError = 7;
        goto endprogram;}
    }
  }
  
  //========================
  //Ouverture du fichier WAV
  //========================
  //Comparaison avec dictionnaire
  if (sDICO)
  {
    save_file = output_file;
    memset (buffer_file, 0x00, 256);
    GetTempPath (256, buffer_file);
    strcat (buffer_file, tmp_file);
    //output_file = (char *)tmp_file;
    output_file = (char *)buffer_file;
    printf ("Read the PCM file and write temporary a MFCC file.\n");     
  }
  
  if (sMFCC || sDICO ||sPlay)
  {
    //Trop peu d'arguments
    if (!input_file || (!output_file && sMFCC)) {
      printf ("Error: Too few argument.\n");
      sError = 1;
      goto endprogram;}
    
    //Trop d'arguments  
    if (sDTW) {
      printf ("Error: Too lot arguments.\n");
      sError = 1;
      goto endprogram;}
      
    //Ouverture du fichier WAV
    if ((sError = WavFile_CheckFormat (input_file, &WAV_wSamples))!=-1)
      goto endprogram; 
      
    //Joue le fichier WAV
    if (sPlay)
    {
      if (!PlayWav_Read (input_file)) 
      {
        printf ("Error: Impossible to play the PCM file!\n");
        sError = 6;
        goto endprogram;
      }
    } 
    
    //MFCC
    if (sMFCC || sDICO)
    {
      //================
      //Traitement audio
      //================
      //Allocation mémoire pour traitement
      if (WavFile_AllocateWindows (WINDOW_SAMPLES, &WAV_wSamples)==false) {
        WavFile_DeleteWindows (&WAV_wSamples);
        printf ("Error: No enough memory!");
        sError = 8;
        goto endprogram;}
  
      //Découpe le fichier WAV en feêntre et traite les échantillons
      WavFile_MakeWindows (WAV_FILE, WINDOW_SAMPLES, &WAV_wSamples);
  
      //Application de MFCC
      WavFile_MFCC (WINDOW_SAMPLES, &WAV_wSamples);
  
      //Enregistrement du spectre dans le fichier MFCC
      WavFile_Create_MFCC_File (&WAV_wSamples, output_file);
  
      //Libère allocations
      WavFile_DeleteWindows (&WAV_wSamples);
      
      //Ferme le fichier PCM
      if (WAV_FILE) fclose (WAV_FILE);         
    }      
  }
  //=============================
  //Dictionnaire
  //=============================  
  if (sDICO)
  {
    FILE *src1, *src2;
    FILE *txtfile;
    
    printf ("Compare PCM file with each files in the dictionnary.\n\n");
    
    //Ouverture MFCC fichier WAV
    memset (buffer_file, 0x00, 256);
    GetTempPath (256, buffer_file);
    strcat (buffer_file, tmp_file);
    if (!(src1 = fopen (buffer_file, "rb"))) {
      printf ("Error: Impossible to open tempory file.\n");
      sError = 1;
      goto endprogram;}
    
    //Création des vecteurs du fichier WAV      
    if (!DTW_CreateVector (src1, &MFCC_wVector1)) {
      DTW_DeleteVector (&MFCC_wVector1);
      fclose (src1);
      goto endprogram;}    

    if (!(txtfile = fopen (save_file, "r"))) {
      printf ("Error: Impossible to open dictionnary file!\n");
      sError = 6;
      goto endprogram;}
    
    //Plus petit DTW
    DTW_value = (float)DTW_INFINITY;
    
    //Efface buffer fichier
    memset (buffer_file, 0x00, 256);
    
    //Parcours fichier    
    while (!feof (txtfile))
    {
      fscanf(txtfile, "%s", buffer);
      //if (feof (txtfile)) break;
      
      //Ouverture MFCC fichier du dico
      if (!(src2 = fopen (buffer, "rb"))) {
        printf ("Error: Impossible to open MFCC file.\n");
        DTW_DeleteVector (&MFCC_wVector1);
        fclose (txtfile);
        fclose (src1);
        sError = 1;
        perror (buffer);
        goto endprogram;}
      
      //Création des vecteurs    
      if (!DTW_CreateVector (src2, &MFCC_wVector2)) {
        DTW_DeleteVector (&MFCC_wVector2);
        DTW_DeleteVector (&MFCC_wVector1);
        fclose (txtfile);
        fclose (src2);
        fclose (src1);
        sError = 1;
        goto endprogram;}
        
      printf ("Fichier %s: Resultat DTW: %0.5f\n", buffer, (DTW_tmp = DTW_Compare (&MFCC_wVector1, &MFCC_wVector2)));
      if (DTW_tmp<DTW_value)
      {
        DTW_value = DTW_tmp;
        GetFileName (buffer, buffer_file);    
      }    
    
      //Efface MFCC_Vector2
      DTW_DeleteVector (&MFCC_wVector2);
      fclose (src2);    
    }     
    
    FILE *out = fopen (OUTPUT_FILE, "w");
    if (DTW_value>DTW_MIDDLE) {
      printf ("No file recognized!\n");
      sError = 1;}
    else 
    {
      printf ("File recognized: '%s' - DTW value: %0.5f\n", buffer_file, DTW_value);
      fputs (buffer_file, out);
    }
    fclose (out);  
          
    //Efface vecteurs
    DTW_DeleteVector (&MFCC_wVector1);
    
    //Ferme fichiers
    fclose (txtfile);
    fclose (src1);      
  }  
  
  //=============================
  //Traitement DTW
  //=============================
  else if (sDTW)
  {
    FILE *src1, *src2;
    
    //Pas de spectres MFCC
    if (!input_file || !output_file) {
      printf ("Error: Too few argument.\n");
      sError = 1;
      goto endprogram;}
   
    //Ouverture des spectres
    if (!(src1 = fopen (input_file, "rb")) || !(src2 = fopen (output_file, "rb"))) {
      printf ("Error: Impossible to open MFCC files.\n");
      sError = 1;
      goto endprogram;}
     
    //Création des vecteurs       
    if (!DTW_CreateVector (src1, &MFCC_wVector1)) {
      DTW_DeleteVector (&MFCC_wVector1);
      fclose (src1);
      goto endprogram;}
    
    if (!DTW_CreateVector (src2, &MFCC_wVector2)) {
      DTW_DeleteVector (&MFCC_wVector2);
      DTW_DeleteVector (&MFCC_wVector1);
    
      //Ferme fichiers
      fclose (src1);
      fclose (src2);
      goto endprogram;}
    
    
    //=====================================
    //TRAITEMENT DTW (Dynamic Time Warping)
    //=====================================
    
    printf ("Resultat DTW: %0.5f\n", DTW_Compare (&MFCC_wVector1, &MFCC_wVector2));
    
    //Efface veteurs
    DTW_DeleteVector (&MFCC_wVector2);
    DTW_DeleteVector (&MFCC_wVector1);
    
    //Ferme fichiers
    fclose (src1);
    fclose (src2);        
  }     
  
  //FIN DU PROGRAMME
  endprogram:
  
  //Ferme la librairie FMOD
  PlayWav_Close ();      
                    
  //Attente appuie sur une touche
  if (sError!=-1)
    system ("PAUSE");
  return sError;
}


//*****************************
//Lecture du format fichier WAV
//*****************************
int WavFile_CheckFormat (char *file, s_wSamples *wSamples)
{
  //Ouverture du fichier
  WAV_FILE = fopen (file,"rb");               //Lecture seule (BINAIRE)
  RiffHeader WAV_HEADER;
  
  if (!WAV_FILE) {
    cerr << "Error: Impossible to open the input file.\n";
    return 1;}
  
  //Enregistrement du Header pour lecture
  fread (&WAV_HEADER, sizeof (RiffHeader), 1, WAV_FILE); 
  
  //Fichier incorrect
  if (memcmp (&WAV_HEADER.FichType,"RIFF", 4) 
   || memcmp (&WAV_HEADER.InfoType, "WAVE", 4)
   || memcmp (&WAV_HEADER.Fmt, "fmt ", 4)) {
    cerr << "Error: Input file is not a PCM file.\n";
    return 2;}
    
  //Format incorrect
  //0x01 = WAVE_FORMAT_PCM
  //0x07 = WAVE_FORMAT_MULAW
  if (WAV_HEADER.wFormatTag!=0x01 && WAV_HEADER.wFormatTag!=0x07) {
    cerr << "Error: Incorrect WAV file format.\n";
    return 2;}
  
  //Vérification de la structure du fichier WAV    
  if (WAV_HEADER.wChannels!=1 || WAV_HEADER.dwSamplesPerSec!=22050 || WAV_HEADER.wBitsPerSample!=16) {
    cerr << "Error: Incorrect format!\n";
    cerr << "Correct format is 'Mono 16 bits at 22050 Hz'.\n";
    return 3;} 
    
  //Section data
  fseek (WAV_FILE, WAV_HEADER.FmtLen+20, SEEK_SET);
  fread (&WAV_HEADER.dwIdentif, 4, 1, WAV_FILE);
  
  //Section 'fact'
  if (memcmp (&WAV_HEADER.dwIdentif, "data", 4))
  {
    fseek (WAV_FILE, 8, SEEK_CUR);
    fread (&WAV_HEADER.dwIdentif, 4, 1, WAV_FILE);
    
    //Vérification section data
    if (memcmp (&WAV_HEADER.dwIdentif, "data", 0)) {
      cerr << "Error: Impossible to open the WAV file.\n";
      return 2;}
  }
  
  //Nombre d'é&chantillons
  fread (&WAV_HEADER.dwNberOfData, sizeof (unsigned long), 1, WAV_FILE);
  
  //Format PCM
  if (WAV_HEADER.wFormatTag==0x01)
    WAV_HEADER.dwNberOfData /= sizeof (short);
  
  //Complètes informations structure
  memcpy (&sHeaderInfo, &WAV_HEADER, sizeof (RiffHeader));                                     
  sNumSamples = WAV_HEADER.dwNberOfData;
                                                                                                                                                                                                                                                              
  return -1;
}


//*****************************************
//Allocation des fenêtres et des structures
//*****************************************
bool WavFile_AllocateWindows (unsigned long nSamples, s_wSamples *wSamples)
{
  //Initialise structure
  sWindow = NULL;
  sMelFilter = NULL;
  sMelFrame = NULL;
  siDCT_Frame = NULL;
  
  unsigned long nWindow = (unsigned long)ceil(((float)sNumSamples/(nSamples/2))-1);
  sNumWindow = nWindow;
  
  if (!(sWindow = (s_Window *)calloc (nWindow, sizeof (s_Window)))) return false; //Fenêtres
  if (!(sMelFilter = (s_MelFilter *)calloc (MELS_FILTERS+1, sizeof (s_MelFilter)))) return false; //Filtres de Mels
  if (!(sMelFrame = (s_data *)calloc (sNumWindow, sizeof (s_data)))) return false; //Mel Frames
  if (!(siDCT_Frame = (s_data *)malloc (sNumWindow*sizeof (s_data)))) return false; //iDCT frames
  
  //Initialisation des données séparées en fenêtres
  for (unsigned long i=0; i<nWindow; i++)
  {
    sWindow [i].data = NULL;
    sWindow [i].FFT_Real = NULL;
    sWindow [i].FFT_Imag = NULL;
    
    sMelFrame [i].data = NULL;
    
    siDCT_Frame [i].data = NULL;
  }
  
  //Initialisation des Filtres de Mels
  for (unsigned long i=0; i<MELS_FILTERS+1; i++)
    sMelFilter [i].data = NULL;
  
  
  //Allocation des données séparées en fenêtres
  for (unsigned long i=0; i<nWindow; i++)
  {
    if (!(sWindow [i].data = (float *)malloc (nSamples*sizeof (float)))) return false;
    if (!(sWindow [i].FFT_Real = (float *)malloc (nSamples*sizeof (float)))) return false;
    if (!(sWindow [i].FFT_Imag = (float *)malloc (nSamples*sizeof (float)))) return false;
    
    if (!(sMelFrame [i].data = (float *)malloc ((MELS_FILTERS+1)*sizeof (float)))) return false;
    memset (sMelFrame [i].data, 0, (MELS_FILTERS+1)*sizeof (float));
    
    if (!(siDCT_Frame [i].data = (float *)malloc ((MELS_FILTERS+1)*sizeof (float)))) return false;
  }
  
  //Allocation des Filtres de Mels
  for (unsigned long i=0; i<MELS_FILTERS+1; i++)
  {
    if (!(sMelFilter [i].data = (float *)malloc ((nSamples/2)*sizeof (float)))) return false;
    memset (sMelFilter [i].data, 0, (nSamples/2)*sizeof (float));
  }
  
  return true;
}


//*********************
//Efface fenêtres
//*********************
void WavFile_DeleteWindows (s_wSamples *wSamples)
{
  unsigned long i;
  
  //Parcours fenêtres
  for (i=0; i<sNumWindow; i++)
  {
    //Efface données FFT
    if (sWindow)
    {
      if (sWindow [i].data) free (sWindow [i].data);
      if (sWindow [i].FFT_Real) free (sWindow [i].FFT_Real);
      if (sWindow [i].FFT_Imag) free (sWindow [i].FFT_Imag);
    }
    
    //Efface MelFrames
    if (sMelFrame && sMelFrame [i].data) free (sMelFrame [i].data);
   
    //Efface iDCT_Frames
    if (siDCT_Frame && siDCT_Frame [i].data) free (siDCT_Frame [i].data);
  }
  
  //Parcours filtres de Mels
  if (sMelFilter)
  {
    for (i=0; i<MELS_FILTERS+1; i++)
      //Efface filtre
      if (sMelFilter [i].data) free (sMelFilter [i].data);
  }
  
  //Efface fenêtres
  if (sWindow) free (sWindow);
  
  //Efface MelFilters
  if (sMelFilter) free (sMelFilter);
  
  //Efface MelFrames
  if (sMelFrame) free (sMelFrame);
  
  //Efface iDCT_Frames
  if (siDCT_Frame) free (siDCT_Frame);
  
  //Efface la structure
  if (wSamples) free (wSamples);
}


//*********************
//Création des fenêtres
//*********************
void WavFile_MakeWindows (FILE *wav_file, unsigned long nSamples, s_wSamples *wSamples)
{  
  //Buffer d'échantillons
  short *SamplesIn = (short *)malloc (nSamples*sizeof (short));
  
  //Création de la fenêtre
  unsigned long SamplesRead = 0;  //Samples déjà lus
  long nReadData = 0;             //Données à ajouter dans le tableau
  bool FirstSample = true;        //Premier échantillon?
   
  //Initialise compteur de fenêtres
  unsigned long nWindow = 0;
        
  //Parcours fichier WAV en fenêtre de NumSamples échantillons
  while (SamplesRead<sNumSamples)
  {
    //Nombre d'échantillons à lire dans le fichier WAV
    nReadData = nSamples;
    
    //Recouvrement des fenêtres de 50% à partir de la 2ème fenêtre, soit avancement de NumSamples/2
    SamplesRead += FirstSample?nSamples:nSamples/2;
    
    //Dépassement du fichier WAV
    if (SamplesRead>wSamples->NumSamples) {
      nReadData = nSamples-(SamplesRead-sNumSamples);
      //Vide buffer échantillons
      memset (SamplesIn, 0, nSamples*sizeof (short));}
    
    //Repositionne pointeur fichier WAV de NumSamples/2 échantillons en arrière
    fseek (wav_file, FirstSample?0:-nSamples/2*sizeof (short), SEEK_CUR);
    
    //Lecture des échantillons
    fread (SamplesIn, nReadData*sizeof (short), 1, wav_file);
    //if (feof (wav_file)) break;
    
    //Marqueur fenetre n>0
    FirstSample = false;
    
    //======================================================================
    //PRETRAITEMENT (Amélioration du signal quasi-stationnaire de la parole)
    //======================================================================
    //Filtre passe haut (préaccentuation)
    /* x(i) = 1-azn^-1 */
    
    //Application fenêtre (Blackman-Harris)
    /*0.355768-0.487396*cos(2*M_PI*x/NumSamples)+0.144232*cos(4*M_PI*x/NumSamples)-0.012604*cos(6*M_PI*x/NumSamples)*/
     
    //Parcours tableau d'échantillons 
    for (unsigned long i=0; i<nSamples; i++)
    {
      //Copie directe en flottant
      sWindow [nWindow].data [i] = (float)SamplesIn [i]/0xFFFF;    //MAX_WORD entre 0 et 0.5???
      
      //Filtre passe haut (préaccentuation)
      //sWindow [nWindow].data [i] = sWindow [nWindow].data [i]?1-PREACCENT_FACTOR/sWindow [nWindow].data [i]:0;
      
      //Fenêtrage du signal (Diminue les fronts montants et calibre périodiquement le signal dans la fenêtre) 
      sWindow [nWindow].data [i] = sWindow [nWindow].data [i]*WindowPonder (i, nSamples);
    }
    
    //Fenêtre suivante
    nWindow++;
  }
  
  //Efface buffer d'échantillons
  free (SamplesIn);
}


//*********************
//Application MFCC
//*********************
void WavFile_MFCC (unsigned long nSamples, s_wSamples *wSamples)
{  
  unsigned long i, j, k;
  
  //===========================================
  //Application de FFT (Fast Fourier Transform)
  //===========================================
  //float real1, imag1, phase1;
  
  //Parcours fenêtres
  for (i=0; i<sNumWindow; i++)
  {
    //Application de FFT
    FFTAudio (nSamples, sWindow [i].data, sWindow [i].FFT_Real, sWindow [i].FFT_Imag);
    
    /*for (unsigned long j=0; j<nSamples/2; j++)
    {
      //db
      //phase = 20*log10 (phase);
      
      //Puissance
      //phase = 2*(real*real+imag*imag)/(nSamples*nSamples);
      
      real1 = wSamples->Window [i].FFT_Real [j];
      imag1 = wSamples->Window [i].FFT_Imag [j];
      phase1 = sqrt (real1*real1+imag1*imag1);///nSamples;
      
      printf ("%f\n", phase1);
    }*/
  }
  
  
  //===========================================
  //Filtres de Mels
  //===========================================
  float MelSize, MelFilterSize;
  
  //Fréquence de maximale 11025 Hz en Mel
  MelSize = FrequencyToMel (sHeaderInfo.dwSamplesPerSec/2);
  
  //Taille d'un filtre
  MelFilterSize = 2*(MelSize/(MELS_FILTERS+1+1));  //Il faut compter MELS_FILTERS+1 filtres
    
  //Création des filtres de Mels
  for (i=0; i<MELS_FILTERS+1; i++)
  {
    //Info filtre
    sMelFilter [i].start_point = i*(MelFilterSize/2);
    sMelFilter [i].end_point = (i*(MelFilterSize/2))+MelFilterSize;
  }
  
  
  float a, b;
  unsigned long U32_a, U32_b;
  
  //Remplis filtres
  //Parcours filtres
  for (i=0; i<MELS_FILTERS+1; i++)
  {
    //Point de départ et d'arrivé en fréquence
    a = MelToFrequency (sMelFilter [i].start_point)/(sHeaderInfo.dwSamplesPerSec/2)*nSamples/2;
    b = MelToFrequency (sMelFilter [i].end_point)/(sHeaderInfo.dwSamplesPerSec/2)*nSamples/2;
    
    U32_a = (unsigned long)floor (a);
    U32_b = (unsigned long)floor (b);
    
    //Création des données du filtre
    for (unsigned long j=U32_a; j<U32_b; j++)
      sMelFilter [i].data [j] = WindowPonder (j-U32_a, U32_b-U32_a);
  }
  
  
  //===========================================
  //Coefficients Cepstraux
  //===========================================
  float real, imag, amp;
    
  //Parcours les frames
  for (i=0; i<sNumWindow; i++)
  {
    //Calcul coefficients Cepstraux
    for (j=0; j<MELS_FILTERS+1; j++)
    {
      //Complexe Ci=x+yi sur nSamples/2
 	  //Calcul coefficients depuis le spectre de puissance
 	  //Formule: ... * sqrt(FFT_Real^2+FFT_Imag^2)^2 = FFT_Real^2+FFT_Imag^2
      
      //Parcours la moitié des valeurs obtenues avec FFT (la seconde partie est en quelque sorte le miroir de la premièer partie)
      for (k=0; k<nSamples/2; k++)
      {
        //Partie Réelle et imaginaire de FFTi
        real = sWindow [i].FFT_Real [k];
        imag = sWindow [i].FFT_Imag [k];
        amp = (real*real)+(imag*imag);
        //amp = sqrt (amp)/sqrt ((float)(SAMPLE_RATE/2));
        
        //Application filtre
        sMelFrame [i].data [j] += sMelFilter [j].data [k]*amp;
      }
      
      //Applique logarithme (Echelle en dB???)
      sMelFrame [i].data [j] = log10 (sMelFrame [i].data [j]);
      //printf ("%f\n", MelFrame [i].data [j]);
    }
  }  
  
  
  //===========================================
  //Calcul iDCT pour chaque frame, créer MFCC
  //===========================================   
  //Application de iDCT
  for (i=0; i<sNumWindow; i++)
  {
    //Application iDCT (Inverse Discret Cosinus Transform)
    iDCT (MELS_FILTERS+1, sMelFrame [i].data, siDCT_Frame [i].data);
    
    /*for (j=0; j<MELS_FILTERS+1; j++)
    {
      printf ("%d\n", (short)(siDCT_Frame [i].data [j]*8192));
    }*/
  }
}


//***********************************************
//Création d'un fichier contenant le spectre MFCC
//***********************************************
bool WavFile_Create_MFCC_File (s_wSamples *wSamples, const char *MFCC_File)
{
  unsigned long i, j;
  short out_value;
  
  FILE *output = fopen (MFCC_File, "wb");
  MFCCHeader MFCC_Header;
  
  //Erreur d'ouverture du fichier
  if (!output) {
    printf ("Error: Impossible to create the output MFCC File.\n");
    return false;}
  
  memcpy (&MFCC_Header.FichType, "MFCC", 4);
  MFCC_Header.NumFilter = MELS_FILTERS;
  MFCC_Header.NumWindow = sNumWindow;
  
  //Ecriture Header
  fwrite (&MFCC_Header, sizeof (MFCCHeader), 1, output);
  
  //Ecriture du spectre
  for (i=0; i<sNumWindow; i++)
  {
    for (j=1; j<MELS_FILTERS+1; j++)
    {
      //Conversion flottant en entier à virgule fixe sur 16 bits signed
      out_value = FloatTo_FIXED_POINT (siDCT_Frame [i].data [j]);
      
      //Ecrit valeur
      fwrite (&out_value, 2, 1, output);
    }
  }

  //Ferme fichier
  fclose (output);

  return true;
}


//***********************************************
//Création des vecteurs
//***********************************************
bool DTW_CreateVector (FILE *input, s_wVector *sVector)
{
  unsigned long i;
  unsigned long NumData, NumDataRead, ReadWindow, ReadFilter;
  short dataS16;
  bool ErrMemory = false;
  
  MFCCHeader MFCC_Header;
  fread (&MFCC_Header, sizeof (MFCCHeader), 1, input);
  
  //Format du fichier correct
  if (memcmp (&MFCC_Header.FichType, "MFCC", 4)) {
    printf ("Error: Incorrect MFCC file format.\n");
    return false;}
  
  //Complète informations vecteurs
  sVector->NumWindow = MFCC_Header.NumWindow;
  sVector->NumFilter = MFCC_Header.NumFilter;
  
  //Allocation des vecteurs
  if (!(sVector->Vector = (s_data *)calloc (MFCC_Header.NumWindow, sizeof (s_data)))) {
    ErrMemory = true;  
    goto endDTW;}
        
  //Initialise vecteurs
  for (i=0; i<MFCC_Header.NumWindow; i++)
    sVector->Vector [i].data = NULL;
  
  //Allocation des données vecteurs
  for (i=0; i<MFCC_Header.NumWindow; i++)
  {
    if (!(sVector->Vector [i].data = (float *)calloc (MFCC_Header.NumFilter, sizeof (float)))) {
      ErrMemory = true;
      goto endDTW;}
  }
  
  //Lecture fichier
  NumData = MFCC_Header.NumWindow*MFCC_Header.NumFilter;
  ReadWindow = 0; ReadFilter = 0;
  NumDataRead = 0;
  
  while (NumDataRead<NumData)
  {
    //Lecture d'un S16 FIXED POINT
    fread (&dataS16, sizeof (short), 1, input);
    if (feof (input)) {
      printf ("Error: During create DTW vectors.\n");
      return false;}
    
    //Remplit vecteur    
    sVector->Vector [ReadWindow].data [ReadFilter] = FIXED_POINT_ToFloat (dataS16);
    
    NumDataRead++;
    if (++ReadFilter>=MFCC_Header.NumFilter)
    {
      ReadFilter = 0;
      ReadWindow++;
    }  
  }
  
  //Fin de la création des vecteurs
  endDTW:
  if (ErrMemory) {
    printf ("Error: No enough memory!\n");
    return false;}
      
  return true;   
}


//***********************************************
//Efface des vecteurs
//***********************************************
void DTW_DeleteVector (s_wVector *sVector)
{
  unsigned long i;
  
  if (sVector->Vector)
  {
    for (i=0; i<sVector->NumWindow; i++)
      if (sVector->Vector [i].data) free (sVector->Vector [i].data);
    free (sVector->Vector);
  }
}

//***********************************************
//Calcul la distance euclidienne
//***********************************************
//Formule: sqrt(sum((x-y)^2));
//     _________________
//     |  n
//     | ___           ²
//     | \   | x - y  |
// r = | /__ |  i   i | 
//    \| i=0
//
float DTW_DistEuclidian (float *x, float *y, unsigned long nFilter1, unsigned long nFilter2)
{
  float r=0;
  for (unsigned long i=0; i<nFilter1 && i<nFilter2; i++)
    r += (x [i]-y [i])*(x [i]-y [i]);
  return sqrt (r);
}

//***********************************************
//Application de DTW
//***********************************************
float DTW_Compare (s_wVector *sVector11, s_wVector *sVector21)
{
  unsigned long nWindow1 = sVector11->NumWindow;
  unsigned long nWindow2 = sVector21->NumWindow;
  unsigned long nFilter1 = sVector11->NumFilter;
  unsigned long nFilter2 = sVector21->NumFilter;
  unsigned long i, j;
  float dist, r;
  
  s_wVector *sVector1 = sVector11;
  s_wVector *sVector2 = sVector21;
  s_wVector *tmp_vec;
  if (nWindow2>nWindow1)
  {
    tmp_vec = sVector1;
    sVector1 = sVector2;
    sVector2 = tmp_vec;
    
    nWindow1 = sVector1->NumWindow;
    nWindow2 = sVector2->NumWindow;
    nFilter1 = sVector1->NumFilter;
    nFilter2 = sVector2->NumFilter;        
  }
  
  float **matrix;
  matrix = (float **)malloc ((nWindow1+1)*sizeof (float));
  for (i=0; i<nWindow1+1; i++)
    matrix [i] = (float *)malloc ((nWindow2+1)*sizeof (float));
        
  for (i=0; i<nWindow1+1; i++)
    matrix [i][0] = DTW_INFINITY;
    
  for (j=0; j<nWindow2+1; j++)
    matrix [0][j] = DTW_INFINITY;
  matrix [0][0] = 0;
  
  for (i=1; i<nWindow1+1; i++)
  {
    for (j=1; j<nWindow2+1; j++)
    {
      //Calcul la distance euclidienne entre 
      //l'item i-1 de x et l'item i de y
      dist = DTW_DistEuclidian (sVector1->Vector [i-1].data, sVector2->Vector [j-1].data, nFilter1, nFilter2);
      matrix [i][j] = min (
                        min (
                          matrix [i-1][j]+dist, 
                            matrix [i-1][j-1]+dist+dist),
                        matrix [i][j-1]+dist);
    }    
  }
  
  r = matrix [nWindow1][nWindow2]/(float)(nWindow1+nWindow2);
  
  for (i=0; i<nWindow1+1; i++)
    free (matrix [i]);
  free (matrix);
  return r;                             
}
