#include "nrutil.h"
#include<stdio.h>
#include <malloc.h>
#include <math.h>
#include <string.h>


#define PI (3.141592654)
#define TWOPI (6.28318530717)
#define EFACTOR (0.1)   /* energy scaling factor */
#define TAPER 0


double log10();
double exp();
double log();
double pow();

/* ========== Parameterizes input speech using LPC analysis ========= */
/**/
int SPEECH_TO_LPC(s,P,nOrd,F_LEN,r,a,engy,gain)
float *engy,*gain;
float *r,*a,*s;
int P,F_LEN,nOrd;
{
float sum;
float *ac;
int i,j;
void CORR();
void SCHUR();
void RCTOA();

 ac = vector(0,nOrd+1);
 
 /* Take the P order autocorrelation of the speech signal s[0..LEN-1] */
 /**/
 CORR(s,ac,nOrd,F_LEN);
 *engy = ac[0];
 if (*engy == 0.0)
   return 1;

 /* Compute the reflection coefficients from the autocorrelation */
 /**/
 SCHUR(ac,r,nOrd);

 /* Convert  the reflection coefficients to predictor coeffs */
 /**/
 RCTOA(r,a,P,nOrd); 

 /* Computes the LP gain */
 /**/
 sum=ac[0]; 
 for (i=1;i<=P;++i)
   sum += a[i]*ac[i];

 if (sum <0.0)
   *gain = 1.0;
 else
 *gain = sqrt(sum);
 
 return 0;
}

/*============== Performs autocorrelation ======================== */

void CORR(s,acf,NP,F_LEN)
float s[];
float acf[];
int NP,F_LEN;
{
int k,j;
float sum;

/* calculate the first NP+1 values of the Auto-correlation function */
/**/
for(k=0; k<=NP ;++k)
 {
   sum=0.0;
   for (j=k ; j<F_LEN; ++j)
      sum += s[j]*s[j-k];

   acf[k] = sum;
 }

}


/* ====== Performs Schur recursion to obtain the reflection coeffs ==*/

void SCHUR(acf,r,NP)
float acf[],r[];
int NP;
{
int n,m,i,j;
float *P,*K;

 P = vector(0,NP);
 K = vector(0,NP);

for(i=1;i<=NP-1;++i)
   K[NP+1-i] = acf[i];

for(j=0;j<=NP;++j)
   P[j] = acf[j];

for (n=1;n<=NP;++n)
{

  if (P[0] < fabs(P[1]))
    {
      for (i=n; i<=NP; ++i)
	r[i] =0.0;
      break;        /* EXIT */
    }
  else
    {
     r[n] = fabs(P[1])/P[0];
     if ( P[1] > 0.0) r[n] = -r[n];
     if (n == NP)
        break;           /* EXIT */
     P[0] += P[1]*r[n];

     for(m=1; m<=NP-n;++m)
     {
        P[m] = P[1+m] + r[n]*K[NP-m+1];
        K[NP-m+1]=K[NP-m+1] + r[n]*P[1+m];
     }
    } /* end of else */
 } /* end of outer for loop */

  } /* end of SCHUR */

/*===========Conversion: from reflection to predictor coeffs ============ */

void RCTOA (k,as,NP,nOrd)
float k[],as[];
int NP,nOrd;
{
int i,j;
float **a;


 a = matrix(1,nOrd+2,1,nOrd+2);

 for (i=1; i<=nOrd; i++ ){    /* perform the  recursion */
   a[i][i]=k[i]; 
  for (j=1;j<i; j++)
      a[j][i]= a[j][i-1]+k[i]*a[i-j][i-1];

}

for (j=1;j<=NP;++j)   /* get the predictor coefficients */
 as[j] = a[j][nOrd];


free_matrix(a,1,nOrd+2,1,nOrd+2);

}


/* ================== Pre-emphasis =================== */
void PREMPHASIZE(beta,s,LEN,FRM)
float beta,*s;
int LEN,FRM;
{
int k;
float temp;
static float prev=0.0;

/*The signal s[0..LEN-1] is applied to a first order FIR preemphasis filter*/
/**/
if (LEN != FRM) prev=s[0];
for(k=0; k<LEN;++k){
  temp = s[k];
  s[k] = temp - beta * prev; 
  prev = temp;
}
 
}

/* ============== Hamming Window ==================== */

void HAMMING(s,LEN)
float *s;
int LEN;
{
float cconst;
int i;
  
  cconst=TWOPI/(LEN-1);
  for (i=0;i<LEN;++i)
     s[i] = s[i]*(0.54-0.46*cos(cconst*i));


}

/* ============ Conversion: from LPC to Cepstrum coeffs =========== */
/**/
void LPC2CEPS(n,nLift,a,cepSave)
int n,nLift;
float *cepSave;
float *a;
{
int i,j;
static fram=0;
float acc;
float *cep;
void CEP_LIFTERING();

  cep = vector(1,n);

      for (i=1;i<=n; ++i){
	acc=0.0;
	for (j=1;j<i;++j)
	   acc += (i-j)*a[j]*cep[i-j];
	cep[i] = -(a[i] +  acc/i);
      }
 
 if (nLift !=0) CEP_LIFTERING(n,nLift,cep);     

   for (i=1; i<=n; ++i)
     cepSave[fram*n+i-1] = cep[i];  /* Save the cepstrum coeffs */
     
fram++;
}

/* ============ Bilinear transform on cepstrum coefficients ===== === */
/**/
void CEPS2BILINEAR(n,BLTord,BLTalpha,nLift,a,cepSave)
int n,nLift,BLTord;
float *cepSave,BLTalpha;
float *a;
{
int i,j;
static fram=0;
float acc;
float *cep,*cc,*implse;
void CEP_LIFTERING();
void BLT_CEPS();



  cep   = vector(1,n);
  cc    = vector (0,BLTord);  /* Bilinear-transformed cepstrum coeffs */
  /* implse= vector(0,n);  */

   /* --- Compute the cesptrum coeffs first ---- */
   /**/
      for (i=1;i<=n; ++i){
	acc=0.0;
	for (j=1;j<i;++j)
	   acc += (i-j)*a[j]*cep[i-j];
	cep[i] = -(a[i] +  acc/i);
      }




 /* if (nLift !=0) CEP_LIFTERING(n,nLift,cep);        */

 /* -- Apply the Bilinear Trasform on cepstrum coeffs -- */
 /**/
 BLT_CEPS(n,cep,BLTord,cc,BLTalpha); 

   for (i=0; i<BLTord; ++i)
     cepSave[fram*BLTord+i] = cc[i];  /* Save the warped coeffs */

free_vector(cep,1,n);
free_vector(cc,0,BLTord);     

fram++;
}


/* =============== Weights the ceptrsum (Cepstrum Liftering) ======== */

void CEP_LIFTERING(n,L,cep)
int n,L;
float *cep;
{
float cconst;
int i;

  cconst = PI/(float)L;
  for (i=1;i<=n;++i)
    cep[i] *= (1.0+0.5*L*sin(cconst*i));


}

/* ============== Regression Coefficients computation ============= */

void REGRESSION_COEF(Rwin,n,numFrames,cep,regres)
int Rwin,n,numFrames;
float *cep,*regres;
{
int i,at,cindx,t,j,k;
int findx;      /* forward index */
int bindx;      /* backward index */
float denom,acc;


 /* -- Compute the denominator of regression coefficients -- */
 /**/
 denom=0.0;
 for (i=1; i<=Rwin; ++i)
    denom += (i*i);
 denom *= 2.0;


 /* --MODIFIED BY Aravind Ganapathiraju on 04/03/98-- */
/* -- Take care of front tail , using first order difference ----- */
/**/
at =0;
for (t=1; t<=Rwin; t++)
  {
    for (j=0;j<n;j++) {
      regres[j] += t*(cep[at+n+j] - cep[j]);
    }    
    at = (t)*n;
  }

for (j=0;j<n;j++) {
  regres[j] /= denom;
}
/*           --> n cepstr <---
        --------+---------+----------+----------+----------   cep() array
	         ^         ^          ^          ^      
		 bindx     at         findx      findx
		  (-n)    current      (+n)      (+2n)
		          location
*/

at = (Rwin -1)*n;
 for (t=Rwin; t<=numFrames-Rwin; t++)  /* -- Full regression formula -- */
   {
     for (k=0;k<n;k++)
       {
	 acc=0.0;
	 findx=k; bindx=findx;
	 for (j=1; j<=Rwin;j++){
	   findx += n;  bindx -= n;
	   if (at==12) {	     
	     acc  += j*( cep[at + findx] - cep[k]);  
	   }
	   else {
	     acc  += j*( cep[at + findx] - cep[at+bindx]);  
	   }
	 } 
	 regres[at+k] = acc/denom;
       }
     at = t*n;  /* Update current location */
   }  /*---  end of full regression -- */

/* -- Take care of end tail , using first order difference ----- */
/**/
 
 for (t=numFrames-Rwin; t<numFrames; t++)
   {    
     for (k=0;k<n;k++)
       {
	 acc=0.0;
	 findx=k; bindx=findx;
	 for (j=1; j<=Rwin;j++){
	   findx += n;  bindx -= n;
	   if ((at + findx) >= (numFrames*n)) {	     
	     acc  += j*( cep[(numFrames-1)*n + k] - cep[at + bindx]);  
	   }
	   else {
	     acc  += j*( cep[at + findx] - cep[at+bindx]);  
	   }
	 } 
	 regres[at+k] = acc/denom;	 
       }
     at = (t+1)*n;
   }
}

		     
/* ============== Normalize the log-energy to be in -1.0 .. 1.0 ==== */

void NORMALIZE_ENGY(energy,numframes)
float *energy;
int numframes;
{
float maxengy=-1000000.0;
int i;

 /* Find maximum energy */
 /**/
  for (i=0;i<numframes; ++i)
    if (energy[i] > maxengy)
      maxengy = energy[i];

 /* Normalize by maximum */
 /**/
  for (i=0;i<numframes; ++i)
     energy[i] = (energy[i] - maxengy)*EFACTOR + 1.0;

}

/* ============== Read in input speech ================== */


void READ_SPEECH(F_LEN,FRM,HDR_TYPE,fpin,s,oengy)
int F_LEN,FRM,HDR_TYPE;
FILE *fpin;
float *s;
float *oengy;
{
static int cnt=1;
int FRAME,number,k,j;
short *shPtr,*shBuf,tmp;
float In_Scale=1.0,engy;
float *buffer;
static float *tail;
void SWAP_BYTES();

if (cnt==1){
   FRAME=F_LEN;
   tail = vector(0,F_LEN-FRM+1);
 }
 else FRAME=FRM;

  buffer = vector(0,F_LEN+1);
  shBuf = (short *)malloc(F_LEN * sizeof(short));

 /* -- Read one frame of speech -- */
 /**/
number= fread(shBuf, sizeof(short), FRAME, fpin);

if (HDR_TYPE==1)               /* - Swap bytes if it's a TIMIT file *.adc-- */
  for (k=0;k<number; ++k){
    tmp = shBuf[k]; 	
    SWAP_BYTES(&tmp);
    shBuf[k]=tmp;
  }


 shPtr = shBuf-1;
 for (k=0;k<FRAME; ++k)
   buffer[k] = *++shPtr * In_Scale;
 

 if (FRM == F_LEN)   /* non-overlapping frames */
   { engy=0.0;
     for (k=0;k<F_LEN;++k){
       s[k]=buffer[k];
       engy = engy + s[k]*s[k];
      }
   }
 else {

   j=0; engy=0.0;
   for (k=0;k<F_LEN;++k){
     if (cnt==1){   /* Take care of 1st frame */
       s[k]=buffer[k];
       if (k>=FRM) 
	 tail[k-FRM]=buffer[k]; 
     }
     else {
       if (k<(F_LEN-FRM))
	 s[k] = tail[k];
       else{
	 s[k]   = buffer[j];
         j++; }
     }
   engy = engy + s[k]*s[k];
   }
 }
*oengy = engy;

if (cnt >1 )
  for (j=FRM; j<F_LEN; ++j)
    tail[j-FRM] = s[j];

 cnt = 2;
}



/* ================ Mel-Filter Bank analysis ============= */

void MEL_CEPSTRUM(NP,F_LEN,N,nOrd,nLift,nLpcMel,LPCspec,SRATE,MelHPF,LFCC,s,engy,
		  Melwt,lbin,FBANK,mfcNorm,NMFCC,cepSave)
int N,NP,F_LEN,*lbin,nOrd,LPCspec,nLpcMel,LFCC,nLift,FBANK,NMFCC;
float *Melwt,*s,*engy,*cepSave,SRATE,MelHPF,*mfcNorm;
{
int i,k,j,Nhalf,bin,ze,WithOFFset=0;
static int irsln=2,binStart=2;
static float rsln;
static int fram=0;
float *y,*bins,magn,real,imag,norm,factor,t1,*freq,engy1,gain;
float acc,*r,*a,*logmagn,denom,liftConst;
static float tpwts[512]; /* Tapering weights */
int SPEECH_TO_LPC();
void FREQ_RESPONSE();
void FFT();


 if (nLift!=0) 
   liftConst = PI/(float)nLift;

 factor = PI/(float)nOrd;
 norm   = sqrt(2.0/(float)nOrd);


 y = vector(1,2*N+2);
 logmagn = vector(1,N/2+1);
 bins = vector(1,nOrd+1);

 Nhalf=N/2;
 k=1;
 acc=0.0;
 for (i=1; i<=N; ++i){
   if (i<=F_LEN){ 
     y[k] = s[i-1];
     y[k+1] =0.0;
     acc += (y[k]*y[k]);  /* Energy computation */
   }
   else {         /* zero pad if necessary */
     y[k]=0.0;
     y[k+1]=0.0;
   } 
   k+=2;
 }

 *engy = acc;

            /* ================================================ */
            /* ------------- Use LPC spectrum ----------------- */
            /* ================================================ */

if (LPCspec) {  /* Apply triang filters to LPC spectrum, rather FFT spec */
  

     a = vector(0,nLpcMel+1); r= vector(0,nLpcMel+1);
     ze    = SPEECH_TO_LPC(s,nLpcMel,nLpcMel,F_LEN,r,a,&engy1,&gain);
     gain=1.0; /* does not matter whether the actual gain is used or not */
               /* since since the spectrum is going to be weighted */
     freq = vector (1,N/2+1); 
     FREQ_RESPONSE(a,gain,freq,N,nLpcMel);

     if (LFCC) {
       if (MelHPF>0.0 && fram==0){ 
	 rsln = SRATE/N;          /* FFT resolution */
	 irsln  = (MelHPF/rsln); /* FFT bin corresponding to MelHPF Hz */
       }
       
       if (MelHPF>0.0)
	 binStart=irsln; 
       
       /* -- Compute cepstrum from IDFT of log-magnitude --- */
       /**/
	for (k=binStart; k<=Nhalf;++k) {
	  magn = freq[k];
	  logmagn[k]=log(magn); 
	}
       
	/*--- Compute cepstrum coeffs ------- */
       /**/ 
       for (i=1; i<=NP; ++i) {
	 acc=0.0;
	 for (k=binStart; k<=Nhalf; ++k) 
	   acc += logmagn[k]*cos(k*TWOPI*i/N);
	  
	 cepSave[fram*NP+i-1] = acc/N;
       }
     }
     else {  /* --Compute cepstrum coeffs from Inv Discrete Cosine Tr---- */

       for (k=2;k<=Nhalf;++k) {
	 magn = freq[k]; 
	 bin  =lbin[k];  
	 t1   =Melwt[k]*magn;
	 if (bin>0) bins[bin] += t1;
	 if (bin<nOrd) bins[bin+1] += magn-t1;
       }
     }

}
            /* ================================================ */
else {      /* -------- Use FFT magnitude spectrum ------------ */
            /* ================================================ */
      FFT(y,N,1);
      

      if (MelHPF>0.0 && fram==0){ 
	rsln = SRATE/N;          /* FFT resolution */
	irsln  = (MelHPF/rsln); /* FFT bin corresponding to MelHPF Hz */
      }

      if (MelHPF>0.0)
	binStart=irsln; 
      
      if (LFCC) {/*---- Use IDFT of log-magnitude to get cepstrum coeff--*/

	/* 
	   if (TAPER && fram==0 && MelHPF>0.0 ){
	   denom=pow(irsln*1.0,4.0);
	   for (i=2;i<irsln; ++i)
	   tpwts[i]=pow(i*1.0,4.0)/denom;
	   }
	   if (TAPER && MelHPF>0.0) binStart=2;
	*/

	
	/* -- Compute cesptrum from IDFT of log-magnitude --- */
	/**/ 
	for (k=binStart; k<=Nhalf;++k) {
	  real=y[2*k-1]; imag=y[2*k];
	  magn = (real*real+imag*imag);  
	  /*
	     if (TAPER && k<irsln && MelHPF>0.0) 
	       logmagn[k]=0.5*log(magn)+log(tpwts[k]); 
	     else
	  */
	    logmagn[k]=0.5*log(magn); 
	}

	/*--- Compute cepstrum coeffs from FFT spectrum ------- */
	/**/
	for (i=1; i<=NP; ++i) {
	  acc=0.0;
	  for (k=binStart; k<=Nhalf; ++k) 
	    acc += logmagn[k]*cos(k*TWOPI*i/N);

	  if (nLift!=0) acc *= (1.0 +0.5*nLift*sin(i*liftConst));
	  cepSave[fram*NP+i-1] = acc/N;
	}

      }
      else {   /*----------- Mel-Cepstrum analysis ---------- */


	for (k=binStart; k<=Nhalf;++k) {
	  real=y[2*k-1]; imag=y[2*k];
	  magn = sqrt(real*real+imag*imag);  
	  /* if (WithOFFset) {  /* ---Shift filterbanks by binStart bins 
	    bin  =lbin[k-binStart+2];  
	    t1   =Melwt[k-binStart+2]*magn; 
	  } 
	  */
	  bin  =lbin[k];  
	  t1   =Melwt[k]*magn; 
	  if (bin>0) bins[bin] += t1;
	  if (bin<nOrd) bins[bin+1] += magn-t1;
	}
      }  /* -- end of if LFCC --- */
      

}

if (!LFCC) { 

        
       for (i=1;i<=nOrd;++i){  
	 acc=bins[i]; 
	 if (acc <1.0) 
	   bins[i]=0.0; 
	   else {
	     if (FBANK || NMFCC)  
	       bins[i]=log(acc/mfcNorm[i]); 
	     else
	       bins[i] = log(acc);
	   }
       }

       if (FBANK) {  /* Dont take the DCT */ 
	 for (i=1; i<=NP; ++i)
	   cepSave[fram*NP+i-1] = bins[i]; /*Save the FBANK coeffs */

       }
       else {

	 for (i=1; i<=NP; ++i) {
	   acc=0.0; t1=i*factor;
	   for (j=1; j<=nOrd;++j)
	     acc += bins[j]*cos(t1*(j-0.5));

	   
	   if (nLift!=0) acc *= (1.0 +0.5*nLift*sin(i*liftConst));
	   cepSave[fram*NP+i-1] = acc*norm;/*Save the mel-cepstrum coeffs */
	 }
       }
}

 free_vector(y,1,2*N+2);
 free_vector(logmagn,1,N/2+1);
 free_vector(bins,1,nOrd+1);

fram++;

}

/* ============= Computes the mel-filter bank weights ================ */

void COMP_FBANK_WEIGHTS(NP,len,SRATE,lbin,wts)
int NP,len,*lbin;
float *wts,SRATE;
{
float *mel,*cenf;
int k;
int b;
void INIT_MEL_SCALE();
void CEN_FREQ();

 mel = vector(1,len+1);
 cenf = vector(1,NP+1);
 
 INIT_MEL_SCALE(len,SRATE,mel);
 CEN_FREQ(NP,mel,len,cenf,lbin);

 for (k=1;k <=len;++k){
   b = lbin[k];
   if (b>0)
       wts[k] = (cenf[b+1]-mel[k])/(cenf[b+1]-cenf[b]);
   else
       wts[k] = (cenf[1] - mel[k])/(cenf[1]);
 /*  printf("%f %f\n",k*0.5*SRATE/len,wts[k]); */
  }

}

/* ========== Initialize the mel-scale =============== */

void INIT_MEL_SCALE(n,SRATE,mel)
int n;
float *mel,SRATE;
{
int k;
float resln,fftsize;

  fftsize = 2.0*n;
  resln = SRATE/(fftsize*700.0); /* Frequency resolution */

  for (k=1;k<=n+1;++k)
    mel[k] = 1127.0*log(1.0+(k-1)*resln);

}

/* =========== Computes the center frequencies ======================= */

void  CEN_FREQ(NP,mel,len,cenf,lowbin)
int NP,len,*lowbin;
float *cenf,*mel;
{
int k,i,bin;
float ms,f;

 ms = mel[len+1];
 for (i=1;i<=NP+1;++i)
   cenf[i] = i*1.0/((float)(NP+1)) * ms;



 bin=1;
 for (i=1; i<=len;++i) {
   f = mel[i];
   while (cenf[bin] < f)
     bin++;
   lowbin[i]=bin-1;
 }

}
/* =============== File handling routines ============================ */

void WRITE_HEADER(num_frames,update,LPC,RegrCoef,Mel,Cepstrum,addEngy,dEngy,NEW,BICEPS,Octave,tparams,fpout)
int num_frames,RegrCoef,Mel,Cepstrum,addEngy,dEngy,
  NEW,LPC,BICEPS,Octave;
short tparams;
double update;    /* Frame update - msecs are implied but are not included */
                  /* in the update number, so update=10.0 means 10.0 msecs */
FILE *fpout;
{
int  SAMPLE_PERIOD;          /* The frame update in 100 nsec units */
double nsec_units=1.0e-7;    /* 100 nsecs */
short ptype=1;               /* Parameterization type */

SAMPLE_PERIOD=update*1.0e-3/nsec_units;

tparams *= 4;   /* Total number of words per feature vector */

/*- Determine type of output parameterization to insert in output header */
/**/
if (LPC && !RegrCoef && !addEngy)
  ptype=1;
if ((Mel || BICEPS || Octave) && !RegrCoef && !addEngy)
  ptype=6;
if (Cepstrum && !RegrCoef && !addEngy)
  ptype=3;

if ((Mel || BICEPS || Octave) && RegrCoef && addEngy)
  ptype=454;
if (Cepstrum && RegrCoef && addEngy)
  ptype=451;

if ((Mel || BICEPS || Octave) && RegrCoef && addEngy && dEngy)
  ptype=0x146;
if (Cepstrum && RegrCoef && addEngy && dEngy)
  ptype=0x143;

if (LPC && (addEngy || dEngy)  && !RegrCoef)
  ptype=65;
if (Cepstrum && ( addEngy || dEngy) && !RegrCoef)
  ptype=67;
if ((Mel || BICEPS || Octave) && (addEngy || dEngy) && !RegrCoef)
  ptype=70;
if ((Mel || BICEPS || Octave) && RegrCoef && !(addEngy || dEngy))
  ptype=0x106;

if (NEW)
  ptype=67;


fwrite(&num_frames,sizeof(int),1,fpout);
fwrite(&SAMPLE_PERIOD,sizeof(int),1,fpout);
fwrite(&tparams,sizeof(short),1,fpout);
fwrite(&ptype,sizeof(short),1,fpout);

}

/* ============ Save all the coefficients in output file ======== */

void  SAVE_COEFFS(num_frames,NP,regrCoef,addEngy,dEngy,cepSave,regres,regres1,
		  Energy,dEnergy,ddEnergy,fpout)
int num_frames,NP,regrCoef,addEngy,dEngy;
float *cepSave,*regres,*regres1,*Energy,*dEnergy,*ddEnergy;
FILE *fpout;
{
float *work1,*work2,*work3,eng,deng,ddeng;
int i;
int k;

work1 = vector(0,NP+1);   /* Scratch, working arrays */
work2 = vector(0,NP+1);
work3 = vector(0,NP+1);

 if (regrCoef) {

   for (k=0;k <num_frames; ++k) {
     for (i=0;i<NP;++i){
       work1[i] =cepSave[k*NP+i];
       work2[i] =regres[k*NP+i];
       work3[i] =regres1[k*NP+i];
     }
     if (addEngy) eng=Energy[k]; 
     if (dEngy) {
       deng=dEnergy[k];
       ddeng=ddEnergy[k];
     }
       
     fwrite(work1,sizeof(float),NP,fpout);
     if (addEngy) fwrite(&eng,sizeof(float),1,fpout); 

     fwrite(work2,sizeof(float),NP,fpout); 
     if (dEngy) {
       fwrite(&deng,sizeof(float),1,fpout); 
     }
     
     fwrite(work3,sizeof(float),NP,fpout); 
     if (dEngy) {
       fwrite(&ddeng,sizeof(float),1,fpout); 
     }
   }
 }
 else
   {

     for (k=0;k <num_frames; ++k) {
       for (i=0;i<NP;++i)
	 work1[i] =cepSave[k*NP+i];
       
       eng=Energy[k];
       
       fwrite(work1,sizeof(float),NP,fpout);
       fwrite(&eng,sizeof(float),1,fpout); 
     }
 }     
     
}


int GetNISTheader (fpin,srate)
FILE *fpin;
float *srate;
{
int num_samples,i,done=0;
char nst[40];


 fscanf(fpin,"%s",nst);   /* Check if it is a valid NIST file */
 if (strcmp(nst,"NIST_1A") == 0) {  /* proceed */
 

   done=0;
   while (!done) {

     fscanf(fpin,"%s",nst);  /* Get next sring */
     if ( strcmp(nst,"sample_count") == 0) {
            fscanf(fpin,"%s",nst);  /* Get the -i string */
	    fscanf(fpin,"%s",nst);  /* Get number of samples */
	    num_samples=atoi(nst);

	  }
     if ( strcmp(nst,"end_head") == 0) {
       printf("ERROR! Could not find NIST keywords in file..\n");
       exit(0);
     }
     if ( strcmp(nst,"sample_rate") == 0) {
       fscanf(fpin,"%s",nst); /* Skip the -i option */
       fscanf(fpin,"%s",nst); /* Get sampling frequency */
       *srate=(float) atoi(nst);  /* Get the Sampling frequency */
       done=1;  /* exit */
     }
   }

  fseek(fpin,1024,0);  /* Skip the 1024 byte header */
 }
else {
  printf("ERROR! Invalid numbers in NIST headers: %s\n",nst);
  exit(0);
}
return(num_samples);
}


/* -------------- NIST header read routine ------------ */

int GetNISTheader_OLD (fpin,srate)
FILE *fpin;
float *srate;
{
int num_samples,i;
char nst[40];


/* --Skip the beginning of header until 'sample_count -i num_samples' --*/
/**/
for (i=1;i<=17;++i)
  fscanf(fpin,"%s",nst);

num_samples=atoi(nst);

for (i=1;i<=3;++i)
 fscanf(fpin,"%s",nst);

*srate=(float) atoi(nst);  /* Get the Sampling frequency */


if (num_samples <160 || num_samples>500000){
  printf("ERROR ! Bad NIST header in num of samples:%d\n\n",num_samples);
  exit(1); }

fseek(fpin,1024,0);  /* Skip the 1024 byte header */

return(num_samples);
}

/* -------------- Byte swapping subroutines ----------------- */
/**/
void SWAP_BYTES(shrt)
short *shrt;
{
   char tmp,*byte;
   
   byte = (char*) shrt;
   tmp = *byte; /* swap MSB with LSB */ 
   *byte = *(byte+1); 
   *(byte+1) = tmp;
}

void SWAP_LONG(word)
int *word;
{
   char tmp,*p;
   
   p = (char*) word;
   tmp = *p; *p = *(p+3); *(p+3) = tmp;
   tmp = *(p+1); *(p+1) = *(p+2); *(p+2) = tmp;
}
/* ---------------- Read label files ------------------------ */
int READ_LBL_FILE(str)
char *str;
{
FILE *fpout;
char str1[40],pl[40],str2[40],str3[40],str4[40];
int beg,end,beg1,end1;


   strcpy(str1,"      ");
   strncpy(str1,str,6); 
   strcat(str1,".phn"); 
   fpout = fopen(str1,"r");
   if (fpout == NULL) { printf("\nERROR! No label file '%s' found\n",str1);
			exit(1);
			}
   else {
     fscanf(fpout,"%d %d %s",&beg,&end,str2);
     fscanf(fpout,"%d %d %s",&beg,&end,str3);
     fscanf(fpout,"%d %d %s",&beg1,&end1,str3);
   }

 if ((end1 - beg-45) <0) 
   printf("ERROR! Not enough frames in: %s\n",str);

/* printf("%d %d: %d %d\n",beg,end,beg+MXFRAMES,end1-(beg+MXFRAMES)); */

 fclose(fpout);


return(beg);

}
/* ------------------- OCTAVE-BANK ANALYSIS ---------------------- */
void OCTAVE_BANK(s,nFilt,NP,N,F_LEN,SRATE,alpha,F0,engy,cepSave)
float *s,F0,*cepSave,SRATE,alpha,*engy;
int NP,nFilt,N,F_LEN;
{
float BW, istart, F1, C, iend, a,denom,cenf,f_engy;
int i,k,Nhalf,vb=0,bin,j;
float rsln,*efreq,real,imag,acc=0.0,*y,*logfilt,*fbank,x,t1;
static int *lbin,fram=0;
static float factor,norm;
void FFT();
 

 Nhalf=N/2;
 
 if (fram==0) {

       factor = PI/(float)nFilt;
       norm   = sqrt(2.0/(float)nFilt);

       efreq = vector(0,nFilt+1);
       lbin  = ivector(0,Nhalf+1);

       a=1.0; denom=1;
       for (i=1; i<nFilt; ++i) {
	a *= alpha;
	denom += a;
       }

       BW = 0.5*SRATE-F0;
       rsln=SRATE/N;  /* Frequency resolution */
       C = BW/denom;  /* Bandwidth of filter 1 */

       /* ---- Determine end frequencies in each band -------- */
       /**/
       istart=F0; F1=C; iend=F0;
       for (i=1; i<=nFilt; ++i) {
	 iend += F1;
	 cenf = istart+0.5*(iend-istart);
	 efreq[i]=iend;     
	 if (vb) 
	   printf("Q[%2d] = %8.3f - %8.3f Hz (%8.3f)\n",i,istart,iend,cenf);
	 istart =iend;
	 F1    *= alpha;
       }

       /* ---- Assign frequency bins to filterbanks --- */
       i=1;
       for (k=1;k<=nFilt; ++k) {
	 while (i*rsln<efreq[k] && i<Nhalf){
	   i++;
	   if (i*rsln>=F0)
	     lbin[i]=k;
	 }
	 if (vb) 
	   printf("i=%3d k=%2d %8.3f\n",i-1,k,efreq[k]);
       }
}

  y = vector(1,2*N+2);
  logfilt = vector(1,nFilt+1);
  fbank   = vector(1,nFilt+1);

 
  k=1;
  acc=0.0;
  for (i=1; i<=N; ++i){
    if (i<=F_LEN){ 
      y[k] = s[i-1];
      y[k+1] =0.0;
      acc += (y[k]*y[k]);  /* Energy computation */
    }
    else {         /* zero pad if necessary */
     y[k]=0.0;
     y[k+1]=0.0;
   } 
   k+=2;
 }

 *engy = acc;

  FFT(y,N,1);

 /* ---- Compute energy in each filterbank ---- */
 /**/
  for (i=1; i<=nFilt;++i) fbank[i]=0.0;

   for (i=1;i<Nhalf ; ++i) {
     bin = lbin[i]; 
     real=y[2*i-1]; imag=y[2*i];
     f_engy = real*real+imag*imag;  
     if (bin>0 && bin<=nFilt)
       fbank[bin] += f_engy;
   }


 /* ---- Take log of filterbank energies -------- */
 /**/
  for (i=1; i<=nFilt; ++i){
    x=fbank[i];
    if (x<1.0)
      x=1.0;
    logfilt[i] = log(x);
  }

 /* -- Take Cosine Transform of log-fbank energies ---- */
 /**/
 for (i=1; i<=NP; ++i) {
  
   acc=0.0; t1=i*factor;
   for (j=1; j<=nFilt;++j)
      acc += logfilt[j]*cos(t1*(j-0.5));
  
   cepSave[fram*NP+i-1] = acc*norm;  /* Save the octave-cepstrum coeffs */
 }

  free_vector(y,1,2*N+2);
  free_vector(logfilt,1,nFilt+1);
  free_vector(fbank,1,nFilt+1);     

fram++;    
  

}