#include <stdio.h>
#include <math.h>

/* assumption: both pattern and template have same number of coefficients */

/************************************************************************/
/* matrixread - allows MFCC data area to be read like a matrix          */
/* Passed:      x: x co-ordinate required (frame)                       */
/*              y: y co-ordinate required (coefficient)                 */
/*              coeffs: number of coefficients per frame                */
/*              *data: pointer to data area                             */
/* Passes:      data value at desired co-ordinate                       */
/************************************************************************/

float matrixread(int x,int y,int coeffs,float *data)
	{
	return(data[(x*coeffs)+y]);
	}


/************************************************************************/
/* ReadMFCC - reads MFCC file into memory                               */
/* Passed:    *inputfile: pointer to opened MFCC file                   */
/*            *coeffs: number of co-efficients in MFCC                  */
/*            *frames: number of frames in MFCC                         */
/* Returns: pointer to (float)allocated memory, or 0 on failure         */
/************************************************************************/

float *ReadMFCC(FILE *inputfile,int *frames,int *coeffs)
	{
	typedef struct 
		{              /* HTK File Header */
		long nSamples;
		long sampPeriod;
		short sampSize;
		short sampKind;
		} HTKhdr;
	HTKhdr hdr;
	float *data;
	int n;

	fread(&hdr, 1, sizeof(HTKhdr), inputfile); 	
		if(hdr.sampKind!= 70)
			{
			fprintf(stderr,"Error: File format not MFCC_E\n");
			return(0);
			}
	data = (float *)malloc(hdr.sampSize*hdr.nSamples);
	n = fread(data,4,(hdr.sampSize/4)*hdr.nSamples,inputfile);
	if(n<(hdr.sampSize/4)*hdr.nSamples)
		{
		fprintf(stderr,"Error: could not read all MFCC values\n");
		free(data);
		return(0);
		}
	*coeffs = hdr.sampSize/4;
	*frames = hdr.nSamples;
	return(data);
	}


/************************************************************************/
/* d - calculates the local distance between pattern frame x and        */
/*     template frame y. Euclidean metric                               */ 
/* Passed:  x: pattern frame index                                      */
/*          coeffs: number of co-efficients in MFCC                     */
/*          *templatedata: pointer to template MFCC data                */
/*          *patterndata: pointer to pattern MFCC data                  */
/* Returns: local distance at desired coordinates in time-time matrix   */
/************************************************************************/

float d(int x,int y, int coeffs, float *templatedata, float *patterndata)
{
  float coeffcost,localcost;
  int a;

  coeffcost=0.0; localcost=0.0;

  for(a=0;a<coeffs;a++) 
    {
    coeffcost=templatedata[(y*coeffs)+a] - patterndata[(x*coeffs)+a];
    localcost=localcost+(coeffcost*coeffcost);
    }
  return(sqrt(localcost));
}


/************************************************************************/
/* min - returns the minimum of three floating point numbers            */
/* Passed:  var1, var2, var3: the three floating point numbers          */
/* Returns: the minimum of the three floating point numbers             */
/************************************************************************/


float min(float var1, float var2, float var3)
{
  if (var1<var2) 
    {
    if (var1<var3) 
      return(var1);
    else 
      return(var3); 
    }
  else 
    {
    if (var2<var3) 
      return(var2); 
    else 
      return(var3); 
    }
}


/************************************************************************/
/* main - speaker independent tests (excluding speaker's data from      */
/*        templates during the recognition process)                     */
/* Algorithm: (see write up for full explanation)                       */
/* Cycles through all template files, grouped by word.                  */
/*   For each of these 'patterns', the recogniser outputs to the screen */
/*   its word estimate using the template files as normal except that   */
/*   templates spoken by the pattern speaker are logically removed from */
/*   the template set.                                                  */
/************************************************************************/

void main(int argc , char **argv)
	{
        /* Begin variable declarations */

        FILE *templateinputfile;    /* file pointer for the template file */
        FILE *patterninputfile;     /* file pointer for the pattern file  */
        int templateframes,patternframes,coeffs;  /* number of template and pattern frames and coefficients */

        float *curCol;              /* current time-time matrix column data     */
       	float *predCol;             /* predecessor time-time matrix column data */

        float *tdata;               /* data for the template file    */
        float *pdata;               /* data for the pattern file     */
        float lowestGlobal;         /* lowest global distance so far */

        int p,i,j;                  /* loop counters */
        int speaker,set,wordnum,pwordnum,pset,pspeaker;  /* loop counters. p indicates pattern */
        int wordindex;              /* word array index for current word estimate */

        char ext[6];                /* file extension .mfcc            */
        char dir[35];               /* directory of templates          */   
        char speakernum[4][2];      /* array holding speaker numbers   */
        char word[13][7];           /* array holding the words spoken  */
        char setnum[8][2];          /* array holding set numbers       */
        char filename[55];          /* final template filename         */
        char pfilename[55];         /* final pattern filename          */

        /* Begin filename construction array definitions */
        strcpy(ext,".mfcc");
        strcpy(dir,"/home/u5snw/com326/templates/");        
        strcpy(speakernum[0],"1");
        strcpy(speakernum[1],"2");
        strcpy(speakernum[2],"3");
        strcpy(speakernum[3],"5");        
        strcpy(word[0],"yes");
        strcpy(word[1],"no");
        strcpy(word[2],"across");
        strcpy(word[3],"down");
        strcpy(word[4],"one");
        strcpy(word[5],"two");
        strcpy(word[6],"three");
        strcpy(word[7],"four");
        strcpy(word[8],"five");
        strcpy(word[9],"six");
        strcpy(word[10],"seven");
        strcpy(word[11],"eight");
        strcpy(word[12],"nine");        
        strcpy(setnum[0],"1");
        strcpy(setnum[1],"2");
        strcpy(setnum[2],"3");
        strcpy(setnum[3],"4");
        strcpy(setnum[4],"5"); 
        strcpy(setnum[5],"6");
        strcpy(setnum[6],"7");
        strcpy(setnum[7],"8");
        
        wordindex=13;    /* default recognised word index. ie Error value */


        printf("Symmetric Dynamic Time Warping Recogniser - Speaker Independent (Def. 2)\n\n");

        /* Read in each 'pattern' data. Abort if the file connot be opened. */
        for(pwordnum=0;pwordnum<13;pwordnum++) 
          {
          printf("\n%s: ",word[pwordnum]);           /* output current word */
          for(pspeaker=0;pspeaker<4;pspeaker++)
            {
            for(pset=0;pset<8;pset++)
              {
              /* construct pattern filename */
              strcpy(pfilename, dir);
              strcat(pfilename, speakernum[pspeaker]);
              strcat(pfilename, "_");
              strcat(pfilename, word[pwordnum]);
              strcat(pfilename, "_");
              strcat(pfilename, setnum[pset]);
              strcat(pfilename, ext);

	      patterninputfile = fopen(pfilename,"r");   /* attempt to open the pattern file */

              /* Only compare if the pattern file could be opened. */
      	      if (patterninputfile!=NULL)
                {
                /* read in the pattern data */
                pdata = ReadMFCC(patterninputfile,&patternframes,&coeffs);

                /* Compare the pattern with each template file. */
                wordindex=13;
                for(speaker=0;speaker<4;speaker++)
                  {
                  if (pspeaker!=speaker)   /* only compare if pattern and template are from different speakers */
                    {
                    for(wordnum=0;wordnum<13;wordnum++) 
                      {
                      for(set=0;set<8;set++)
                        {
                        /* construct template filename */
                        strcpy(filename, dir);
                        strcat(filename, speakernum[speaker]);
                        strcat(filename, "_");
                        strcat(filename, word[wordnum]);
                        strcat(filename, "_");
                        strcat(filename, setnum[set]);
                        strcat(filename, ext);
  
                        if (strcmp(pfilename,filename)!=0)   /* don't compare if files are same */
                          {
         	          templateinputfile = fopen(filename,"r");  /* attempt to open template file */

                          /* Only compare if the template could be opened. */
                          if (templateinputfile!=NULL)
                            {
                            /* read in template data */
                            tdata = ReadMFCC(templateinputfile,&templateframes,&coeffs);

                            /* allocate memory for the 1st two columns of the time-time matrix */
                            predCol = (float *)malloc(sizeof(float)*templateframes);
                            curCol  = (float *)malloc(sizeof(float)*templateframes);

                            /* calculate 1st column of time-time matrix */
                            predCol[0]=d(0,0,coeffs,tdata,pdata);
                            for(p=1;p<templateframes;p++) 
                              predCol[p]=predCol[p-1]+d(0,p,coeffs,tdata,pdata); 
  
                            /* calculate all subsequent columns */
                            for(i=1;i<patternframes;i++) 
                              { 
                              curCol[0]=predCol[0]+d(i,0,coeffs,tdata,pdata);  
                              for(j=1;j<templateframes;j++)  
                                curCol[j]=(min(predCol[j],predCol[j-1],curCol[j-1])+d(i,j,coeffs,tdata,pdata)); 
                              free(predCol);  
                              predCol=curCol;
                              curCol = (float *)malloc(sizeof(float)*templateframes); 
                              }

                            /* if this template is a better match than previous comparisons, make it   */
                            /* the new word estimate. If this is the first time round the loop, merely */
                            /* assign the best estimate so far to be this template's word.             */

                            if (wordindex==13) /* ie first time round loop(s) */
                              {
                              lowestGlobal=predCol[templateframes-1];  /* top right of time-time matrix */
                              wordindex=wordnum; 
                              }
                            else 
                              {
                              if (predCol[templateframes-1] < lowestGlobal)
                                {
                                lowestGlobal=predCol[templateframes-1];
                                wordindex=wordnum;
                                }
                              }
                            
                            /* garbage collect the 2 columns of the time-time matrix ready for the next loop, or termination */
                            free(curCol);
                            free(predCol);

                            /* garbage collect the memory used for the templates data */
                            free(tdata);  
                            fclose(templateinputfile);  /* close the template file */
                            } 
                          }
                        }
                      }    
                    }
                  }

                /* garbage collect the pattern's data */
                free(pdata);
                fclose(patterninputfile);   /* close the pattern file */

                /* now output the word estimate of the pattern file */
                if (wordindex!=13) 
                  printf(" %s ",word[wordindex]);
                else
                  printf("Error locating template files in %s\n",dir);
                }
              }
            }
          }
        }