#include "fft.h"

int main(int argc , char *argv[])
{  
	
int j,i,k=0,NUM=0;
double s[26];	
double c[16];
double H[256];




double f[27];
int f_r[27];

double mel;
double f1=1125*log(1+(double)20/(double)700);//求出mel濾波器的25個中心頻率
//double f1=0;
double f2=1125*log(1+(double)5500/(double)700);
for(i=0; i<27; i++)                   
{
mel=(f1+i*(f2-f1)/(double)27)/(double)1125;           
f[i]=((double)256/(double)11025)*700*(exp(mel)-1);   
f_r[i]=round(f[i]); 
}        



  


char fname[100],fname1[100];
char str[256];
FILE *fp;
FILE *fp_out;
FILE *fp_in;







if(!(fp=fopen(argv[1],"r")))
 {
  printf("Can not open file: %s\n",argv[1]);
  exit(0);
 }


strcpy(fname1,"D:\\feature\\sss\\clean10ci");//結果文件存放的路徑
//strcat(fname1,str);
fp_out = fopen(fname1, "w");
while( fscanf(fp, "%s", str ) != EOF)//語音文件循環開始
 {
  strcpy(fname,"D:\\audio\\test7ren\\clean10ci\\");//原始語音存放的路徑
  strcat(fname,str);
  if(!(fp_in=fopen(fname,"r")))
	{
		printf("Can not open file:%s\n",fname);
		exit(0);
	}



double **mfcc;
mfcc=new double*[200];
for( i=0;i<200;i++)
{
 mfcc[i]=new double[256];
}	
double **ci;
ci=new double*[200];
for( i=0;i<200;i++)
{
  ci[i]=new double[iNUM];
}






float tmp;
int iLen;
iLen=0;
double *PsSpeech=new double[MAX_SPEECH_LEN];
double *Speech=new double[MAX_SPEECH_LEN];
while(fscanf(fp_in,"%f",&tmp)!=EOF)
    {
     PsSpeech[iLen]=(double)tmp;
	 iLen++;
	}



printf("NUM=%d\n",NUM);
NUM++;
//printf("iLen=%d\n",iLen);
fclose(fp_in);




for(i=1;i<iLen;i++)  //預加重
{
 Speech[i-1]=PsSpeech[i]-0.95*PsSpeech[i-1];
 }



int len=0;
int frame_num=0;
int Ps_num=0;
while(!((len+256)>(iLen-2)))//分幀循環開始
{
for( i=0;i<256; i++)  
{   
  mfcc[frame_num][i]= Speech[Ps_num]*(0.54-0.46*cos(2*i*PI/(double)255));
  Ps_num++;
}

len=len+128;
Ps_num=len;
frame_num++;
}//分幀循環結束
delete[] PsSpeech;
delete[] Speech;
/*k=0;
for(i=0; i<frame_num; i++)
{
for(j=0; j<256; j++)
{
printf("mfcc=%f\n",mfcc[i][j]);

printf("k=%d\n",k);
k++;
}
}*/





for(i=0; i<frame_num; i++)//幀循環開始
{

double xxi[256];
for(j=0;j<256;j++)
{xxi[j]=0;}

splfft(mfcc[i],xxi,256,-1);//經過FFT變換

/*
for(j=0; j<256; j++)
{
printf("mfcc[%d][%d]=%f\n",i,j,mfcc[i][j]);

}*/






for(j=0; j<26; j++)//將s[]清零,以免上次的值干擾;s[]保存24個mel濾波器的輸出
{
 s[j]=0;
}
                    




for(j=1; j<26; j++)//計算濾波器的值,并將fft的值通過濾波器累加,形成24個s[]值                   
{ 
for (k=f_r[j-1]; k<= f_r[j+1]; k++)                  
{
	if(k<=f_r[j])
	{	H[k]=((double)(k-f_r[j-1]))/((double)(f_r[j]-f_r[j-1]));
	//printf("liang=%d\n",k);
	}
    else 
	{   H[k]=((double)(f_r[j+1]-k))/((double)(f_r[j+1]-f_r[j]));
	//printf("zhou=%d\n",k);
	}

s[j]+=H[k]*(mfcc[i][k]*mfcc[i][k]);

} 

s[j]=log(s[j]);
}






for(j=0; j<16; j++)//對c[n]清零,以消除上次值的影響; c[n]存放的是一幀數據中提取的16個特征
	{
	 c[j]=0;
	}


for(j=0;j<16;j++)
{
for(int k=1; k<26; k++)
{
   //double fg=0;
   c[j]+=s[k]*(cos(j*PI*(k-0.5)/26));
   //c[j]=c[j]+fg;
}
 //  c[j]=(1+6*sin(j*PI/16))*c[j];//對倒譜進行提升
 
   ci[i][j]=sqrt((double)2/(double)26)*c[j];
  
}

/*for(j=22;j<26;j++)
{
for(int k=1;k<26;k++)
{
   double fg=0;
   fg=s[k]*(cos(j*PI*(k+0.5)/26));
   c[j]=c[j]+fg;
}
 //  c[j]=(1+6*sin(j*PI/16))*c[j];//對倒譜進行提升

   ci[i][j-10]=c[j];
  
}*/

}//幀循環的結束
/*
 double DON1=0,DON2=0;//求取動態特征


for(i=0; i<frame_num; i++)
{
 if(i<3)
 {
  for(j=0; j<16; j++)
  {
   ci[i][j+16]=ci[i+1][j]-ci[i][j];
  }
 }
 else if(i>frame_num-4)
 {
  for(j=0; j<16; j++)
  {
   ci[i][j+16]=ci[i-1][j]-ci[i][j];
  }
 }
 else
 {
  for(j=0; j<16; j++)
  {      
      for(k=0; k<3; k++)
	  {
       DON1+=(ci[i+k][j]-ci[i-k][j])*k*(0.54-0.46*cos(2*k*PI/(double)6));
   
	   DON2+=k*k*(0.54-0.46*cos(2*k*PI/(double)6));
	  }
	  ci[i][j+16]=DON1/DON2;
      DON1=0;
	  DON2=0;
  }
 }
}//求動態特征結束

*/

/*
strcpy(fname0,"F:\\mfcc\\feature\\clean\\TT");
strcat(fname0,str);
fp_outl=fopen(fname0, "w");
for(j=0; j<i; j++)
{
for(k=0; k<16; k++)
{
fprintf(fp_outl,"%f\n", ci[j][k] );
}
fprintf(fp_outl,"\n");
}
fclose(fp_outl);
*/









double *cmfcc=new double [1024];
TimeNormalize(ci,  frame_num,  cmfcc, 64);
AmpNormalize( cmfcc,  1024);



fprintf(fp_out,"%d ",((NUM-1)%10+1));//n詞就改成n
for( i=0 ; i<1024 ;i++ )
{
//	if( !(i%16) ) fprintf(fp_out,"\n");
    fprintf(fp_out,"%d%c%f ",i+1,':',cmfcc[i]);
}
fprintf(fp_out,"\n");

for(j=0; j<200; j++)
delete [] ci[j];
delete [] ci;
ci=NULL;
for(j=0; j<200; j++)
delete [] mfcc[j];
delete [] mfcc;
mfcc=NULL;


delete [] cmfcc;
cmfcc=NULL;

}//語音文件循環結束
fclose(fp_out);
fclose(fp);
return 0;
}



int round(double d)
{
 int f;
  f=(int)d;
  if((double)(f+0.5)>d)
        f=f;
   else
    	f=f+1;   
  return f;
}


void bubble(double *items,int count)
{
int a,b;
double t;
for(a=1; a<count; ++a)
for(b=count-1; b>=a; --b)
{
if(items[b-1]>items[b])
{
 t=items[b-1];
 items[b-1]=items[b];
 items[b]=t;
}
}
}