snr.c

我的语音直接波形矢量量化器,一个简单的例子.没有考虑有噪信道下的VQ-er设计

/*************
       SNR algorithm.

       1).Seg_SNR.
       2).Glob_SNR.

    Input file is "*.wav".
	 When not,the first 0x2c byte should be included.

	L.JY   1997.4.
	 2000.8.
					  *******************/

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

#define NFR 256

//  Compute average segmental SNR for segments above threshold
  void seg_snr(char *in1,char *in2,float *seg, short int *nlow,short int *nseg)
     {
	  FILE *fp_i1,*fp_i2;
	  int cnt,cnt1,cnt2;
	  int i;

	  short int ix1[NFR],ix2[NFR];
	  float x[NFR],y[NFR],xp,yp,sum,xth,maxsnr;

  if((fp_i1=fopen(in1,"rb"))==NULL)
   {  printf("Can not open file %s!\n",in1);
      exit(1); }
  if((fp_i2=fopen(in2,"rb"))==NULL)
   {  printf("Can not write file %s!\n",in2);
      exit(1); }

     //When input data file is speech file in format "*.wav".
      fseek(fp_i1,0x2c,SEEK_SET);
      fseek(fp_i2,0x2c,SEEK_SET);

	xth=10000.0;
	maxsnr=200.0;
	sum=0;
	*nlow=0;
	*nseg=0;

    while( (!feof(fp_i1)) && (!feof(fp_i2) ) )
    {
       cnt1=fread(ix1,sizeof(short int),NFR,fp_i1);
       cnt2=fread(ix2,sizeof(short int),NFR,fp_i2);
       if(cnt1<=cnt2)
	cnt=cnt1;
       else
	cnt=cnt2;

       for(i=0;i<cnt;i++)
	{ x[i]=(float)ix1[i];
	  y[i]=(float)ix1[i]-(float)ix2[i]; }

	xp=0;
	yp=0;
	for(i=0;i<cnt;i++)
	 { xp+=x[i]*x[i];
	   yp+=y[i]*y[i];  } //for i
	xp=xp/cnt;
	yp=yp/cnt;

	if(xp<=xth)
	 (*nlow)++;
	else
	 { (*nseg)++;
	   if(yp>0)
	    sum+=10.0*log10(xp/yp);
	   else
	    sum+=maxsnr;
		 } //else

      } //while read in data

	*seg=sum/(*nseg);

	     fclose(fp_i1);
	     fclose(fp_i2);

	 }


//  Compute Global SNR
  void glob_snr(char *in1, char *in2, float *glb)
     {
	  FILE *fp_i1,*fp_i2;
	  int cnt,cnt1,cnt2;
	  int i;

	  short int ix1[NFR],ix2[NFR];
	  float x[NFR],y[NFR],xp,yp;

  if((fp_i1=fopen(in1,"rb"))==NULL)
   {  printf("Can not open file %s!\n",in1);
      exit(1); }
  if((fp_i2=fopen(in2,"rb"))==NULL)
   {  printf("Can not write file %s!\n",in2);
      exit(1); }

      fseek(fp_i1,0x2c,SEEK_SET);
      fseek(fp_i2,0x2c,SEEK_SET);

      xp=0;
      yp=0;

    while( (!feof(fp_i1)) && (!feof(fp_i2) ) )
    {
       cnt1=fread(ix1,sizeof(short int),NFR,fp_i1);
       cnt2=fread(ix2,sizeof(short int),NFR,fp_i2);
       if(cnt1<=cnt2)
	cnt=cnt1;
       else
	cnt=cnt2;

       for(i=0;i<cnt;i++)
	{ x[i]=(float)ix1[i];
	  y[i]=(float)ix1[i]-(float)ix2[i]; }

	for(i=0;i<cnt;i++)
	 { xp+=x[i]*x[i];
	   yp+=y[i]*y[i]; }

      } //while read in data

	if(xp<=0)
	 *glb=-1;
	else
	 if(yp<=0)
	  *glb=-2;
	 else
	  *glb=10.0*log10(xp/yp);

       fclose(fp_i1);
       fclose(fp_i2);

	 }