97wave.c

给予TMS320C6713浮点DSP的97小波语音信号分析算法实现

/**********************************************************
development by [e-mail]
**********************************************************/

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

#define F_EXTPAD 4
#define D_EXTPAD 2

typedef unsigned int       DWORD;
typedef int                 BOOL;
typedef unsigned char       BYTE;
typedef unsigned short      WORD;
typedef float               FLOAT;

/*  float 9/7
you can get the coefficient by using [a,b,c,d]=WFILTERS('bior4.4') in matlab
********************************/
double h[5] ={ 0.85269867900889 ,  0.37740285561283 , -0.11062440441844,  -0.02384946501956 ,  0.03782845550726 };
double g[4] ={ -0.78848561640558 ,  0.41809227322162 , 0.04068941760916,  -0.06453888262870  };
double h_b[4] ={ 0.78848561640558 ,  0.41809227322162 , -0.04068941760916,  -0.06453888262870 };
double g_b[5] ={-0.85269867900889 ,  0.37740285561283 , 0.11062440441844,  -0.02384946501956 ,  -0.03782845550726  };

/*********  float 5/3 **********
you can get the coefficient by using [a,b,c,d]=WFILTERS('bior2.2') in matlab
**********************************
double h[5] ={  1.06066017177982 ,  0.35355339059327  ,-0.17677669529664 ,0};
double g[4] ={ -0.70710678118655 ,  0.35355339059327 , 0 };
double h_b[4] ={ 0.70710678118655 ,  0.35355339059327, 0 };
double g_b[5] ={-1.06066017177982 ,  0.35355339059327,   0.17677669529664,0 };
*/
static float   *x_alloc = NULL;	// **

short *sample_data;

float data_f[128];

int readfile(char *filename,short **speech,int *sample_length)
{
    FILE *fp;
    BYTE id[5];
    DWORD data_id;
    DWORD size; 
    short format_tag, block_align, bits_per_sample,channels;        // 16 bit values
    DWORD format_length,avg_bytes_sec,sample_rate;    // 32 bit values


	if((fp = fopen(filename,"rb"))==NULL)
	{
		printf("can not open file!");
	    return 0;
	}

    fread(id, sizeof(BYTE), 4, fp);         //"RIFF"标志，辨别码
    id[4]='\0';
	if (!strcmp((char*)id, "RIFF"))         
    { 
        fread(&size, sizeof(DWORD), 1, fp); //文件长度
        fread(id, sizeof(BYTE), 4, fp);     //"WAVE"标志，格式辨别码
        id[4]='\0';
	    if (!strcmp((char *)id,"WAVE")) 
        {
            fread(id, sizeof(BYTE), 4, fp);                //"fmt"标志，FORMAT块
            
			fread(&format_length, sizeof(DWORD),1,fp);     // fmt块长度16或20 			  
            fread(&format_tag, sizeof(short), 1, fp);      //编码格式，如WAVE_FORMAT_PCM，WAVE_F0RAM_ADPCM等,01H为PCM文件格式
		    
		    fread(&channels, sizeof(short),1,fp);           //声道数：1单声道，2双声道
			printf("声 道 数 = %d(个)\n",channels);
			
			fread(&sample_rate, sizeof(DWORD), 1, fp);      //采样率8000
			printf("采 样 率 = %d(Hz)\n",sample_rate);
			
			fread(&avg_bytes_sec, sizeof(DWORD), 1, fp);   //波形音频数据传送速率=通道数×每秒数据位数×每样本的数据位数／8
			
			fread(&block_align, sizeof(short), 1, fp);     //一个采样值2BYTE
			fread(&bits_per_sample, sizeof(short), 1, fp); //采样精度（幅度）(16 bit) 
			printf("采样精度 ＝%d(位)\n",bits_per_sample);
			
			fread(&data_id, sizeof(DWORD), 1, fp);                //数据标记符＂data＂,data块
			
			fread(sample_length, sizeof(DWORD), 1, fp);    //语音数据的长度
			printf("数据长度 = %d(字节)\n",*sample_length);

			*sample_length = *sample_length/2;             //speech是SHORT=2BYTE,信号数据共有*sample_length/2个双字节数；
   			*speech = (short *) malloc (sizeof(short) * (*sample_length)); //set aside sound buffer space
			fread(*speech, sizeof(short), (*sample_length), fp); //read in our whole sound data chunk
		}
	}
	
	fclose(fp);
    return 0;
}


unsigned char **Alloc_Char_Matrix(short r, short c)
{
   unsigned char  *x, **y;
   short           n;
   // calloc
   x = (unsigned char *) calloc((long) r * c, sizeof(unsigned char));
   y = (unsigned char **) calloc(r, sizeof(unsigned char *));
   for (n = 0; n <= r - 1; n++)
       y[n] = &x[(long) c * n];
   return y;
}

short **Alloc_Short_Matrix(short r, short c)
{
   short          *x, **y;
   short           n;
   // calloc
   x = (short *) calloc((long) r * c, sizeof(short));
   y = (short **) calloc(r, sizeof(short *));
   for (n = 0; n <= r - 1; n++)
      y[n] = &x[(long) c * n];
   return y;
}

float **Alloc_Float_Matrix(short r, short c)
{
   float          *x, **y;
   short           n;
   // calloc
   x = (float *) calloc((long) r * c, sizeof(float));
   y = (float **) calloc(r, sizeof(float *));
   for (n = 0; n <= r - 1; n++)
       y[n] = &x[(long) c * n];
   return y;
}

void File2_Char_Matrix(unsigned char *f_name, unsigned char **x, short r, short c)
{
   short           i, j;
   FILE           *f_in;
   f_in = fopen(f_name, "rb");
   for (i = 0; i < r; i++)
     for (j = 0; j < c; j++)
	 fread(&x[i][j], sizeof(unsigned char), 1, f_in);
   fclose(f_in);
}

void  Char_Matrix_2File(unsigned char *f_name, unsigned char **x, short r,short c)
{
   short           i, j;
   FILE           *f_out;
   f_out = fopen(f_name, "wb");
   for (i = 0; i < r; i++)
     for (j = 0; j < c; j++)
	 fwrite(&x[i][j], sizeof(unsigned char), 1, f_out);
   fclose(f_out);
}


void File2_Short_Matrix(unsigned char *f_name, short **x, short r, short c)
{
   short           i, j;
   FILE           *f_in;
   f_in = fopen(f_name, "rb");
   for (i = 0; i < r; i++)
     for (j = 0; j < c; j++)
	 fread(&x[i][j], sizeof(short), 1, f_in);
   fclose(f_in);
}

void  Short_Matrix_2File(unsigned char *f_name, short **x, short r,short c)
{
   short           i, j;
   FILE           *f_out;
   f_out = fopen(f_name, "wb");
   for (i = 0; i < r; i++)
      for (j = 0; j < c; j++)
	 fwrite(&x[i][j], sizeof(short), 1, f_out);
   fclose(f_out);
}


void File2_Float_Matrix(unsigned char *f_name, float **x, short r, short c)
{
   short           i, j;
   FILE           *f_in;
   f_in = fopen(f_name, "rb");
   for (i = 0; i < r; i++)
     for (j = 0; j < c; j++)
	 fread(&x[i][j], sizeof(float), 1, f_in);
   fclose(f_in);
}

void  Float_Matrix_2File(unsigned char *f_name, float **x, short r,short c)
{
   short           i, j;
   FILE           *f_out;
   f_out = fopen(f_name, "wb");
   for (i = 0; i < r; i++)
     for (j = 0; j < c; j++)
	 fwrite(&x[i][j], sizeof(float), 1, f_out);
   fclose(f_out);
}



void  forwardf97(float *x_in, short N)
{
   short           i, n, half;
   float          *x, *r, *d;

   x = x_alloc + F_EXTPAD;
   memcpy(x, x_in, sizeof(float) * N);

   for (i = 1; i <= F_EXTPAD; i++)
     {
      x[-i] = x[i];	  
      x[(N - 1) + i] = x[(N - 1) - i];
     }   

   half = N >> 1;

   r = x_in;
   d = x_in + half;

   for (n = half; n--;)
   {

	  *r++ = h[4] * (x[4] + x[-4]) + h[3] * (x[3] + x[-3]) + h[2] * (x[2] + x[-2]) +
		 h[1] * (x[1] + x[-1]) + h[0] * x[0];

	  x++;

	  *d++ = g[3] * (x[3] + x[-3]) + g[2] * (x[2] + x[-2]) + g[1] * (x[1] + x[-1]) +
		 g[0] * x[0];

	  x++;
   }

}

void  inversef97(float *x, short N)
{
   short           i, n, half;
   float          *r, *d;

   half = N / 2;

   r = x_alloc + D_EXTPAD;
   d = x_alloc + D_EXTPAD + half + D_EXTPAD + D_EXTPAD;
   memcpy(r, x, half * sizeof(float));
   memcpy(d, x + half, half * sizeof(float));

   for (i = 1; i <= D_EXTPAD; i++)
   {
	  r[-i] = r[i];
	  r[(half - 1) + i] = r[half - i];
	  d[-i] = d[i - 1];
	  d[(half - 1) + i] = d[(half - 1) - i];
   }

   for (n = half; n--;)
   {

	  *x++ = h_b[0] * r[0] + h_b[2] * (r[1] + r[-1]) + g_b[3] * (d[1] + d[-2]) +
		 g_b[1] * (d[0] + d[-1]);

	  *x++ = h_b[1] * (r[1] + r[0]) + h_b[3] * (r[2] + r[-1]) + g_b[4] * (d[2] + d[-2]) +
		 g_b[2] * (d[1] + d[-1]) + g_b[0] * d[0];

	  d++;
	  r++;
   }

}

void  f97_1D(float *rows, short width, short levels, short inverse)
{
   short           x, y, w, h, l;
   long            yw;
   float         *buffer, *rows_ptr;

   /*Check the dimensions for compatability.   */

   if (width % (1 << levels))
   {
	  printf("width and height must be divisible by 2^levels");
	  exit(10);
   }

   if ((x_alloc =malloc(sizeof(float) * (width + F_EXTPAD + F_EXTPAD))) == NULL)
   {
	  printf("malloc failed");
	  exit(10);
   }


   /*Compute the rows wavelet transform.  */

   if (!inverse)
   {	/*forward transform. */

	  for (l = 0; l < levels; l++)
	  {
		 w = width >> l;
		 /*Rows. */	
		 forwardf97(rows, w);		
	  }

   } else
   {

	  for (l = levels - 1; l >= 0; l--)
	  {
		 w = width >> l;
		 /* Columns. */
		 inversef97(rows, w);
	  }
   }

   free(x_alloc);
   x_alloc = NULL;
}

void  f97_2D(float **rows, short width, short height, short levels, short inverse)
{
   short           x, y, w, h, l;
   long            yw;
   float         *buffer, *rows_ptr;

   /*Check the dimensions for compatability.   */

   if (width % (1 << levels) || height % (1 << levels))
   {
	  printf("width and height must be divisible by 2^levels");
	  exit(10);
   }

   if ((x_alloc =malloc(sizeof(float) * (width + height + F_EXTPAD + F_EXTPAD))) == NULL)
   {
	  printf("malloc failed");
	  exit(10);
   }

   /* Allocate a work array (for transposing columns) */

   if ((buffer = malloc(sizeof(float) * height)) == NULL)
   {
	  printf("malloc failed");
	  exit(10);
   }

   /*Compute the rows wavelet transform.  */

   if (!inverse)
   {	/*forward transform. */

	  for (l = 0; l < levels; l++)
	  {
		 w = width >> l;
		 h = height >> l;

		 /*Rows. */

		 for (y = 0; y < h; y++)
			forwardf97(rows[y], w);

		 /* Columns. */

		 rows_ptr = rows[0];
		 for (x = 0; x < w; x++)
		 {
			for (y = 0, yw = 0; y < h; y++, yw += width)
			   buffer[y] = rows_ptr[yw];
			forwardf97(buffer, h);
			for (y = 0, yw = 0; y < h; y++, yw += width)
			   rows_ptr[yw] = buffer[y];
			rows_ptr++;
		 }
	  }

   } else
   {

	  for (l = levels - 1; l >= 0; l--)
	  {
		 w = width >> l;
		 h = height >> l;

		 /* Columns. */

		 rows_ptr = rows[0];
		 for (x = 0; x < w; x++)
		 {
			for (y = 0, yw = 0; y < h; y++, yw += width)
			   buffer[y] = rows_ptr[yw];
			inversef97(buffer, h);
			for (y = 0, yw = 0; y < h; y++, yw += width)
			   rows_ptr[yw] = buffer[y];
			rows_ptr++;
		 }

		 /* Rows. */

		 for (y = 0; y < h; y++)
			inversef97(rows[y], w);
	  }
   }

   free(x_alloc);
   x_alloc = NULL;
   free(buffer);
}

float PSNR(unsigned char *in_file,unsigned char *out_file,short x_size,short y_size)
 {
  short i,j;
  float MSE=0.0;
  unsigned char **in_img,**out_img;
  in_img = Alloc_Char_Matrix(x_size, y_size);
  out_img = Alloc_Char_Matrix(x_size, y_size);

  File2_Char_Matrix(in_file, in_img, x_size, y_size);
  File2_Char_Matrix(out_file,out_img, x_size, y_size);

  for(i=0;i<x_size;i++)
     for(j=0;j<y_size;j++)
      MSE+=((float)in_img[i][j]-out_img[i][j])*(in_img[i][j]-out_img[i][j]);

   return(10*log10(255.0*255*x_size*y_size/MSE));
 }


short main()
{
	short           i,j;
   float           x[8]={1.5,25.0,693.0,0.58,1.35,15.2,2.35,63.2};
   int             datalen;
   /*
   char            in_f[60], out1_f[60], out2_f[60];
   unsigned char   **char_img;
   float           **in_img;
   short           **out_img;
   short           i,j,x_size, y_size, levels;


   strcpy(in_f, "lena512.raw");
   strcpy(out1_f, "data11.raw");
   strcpy(out2_f, "out.raw");
   x_size = 8;
   y_size = 8;
   levels = 2;


   in_img = Alloc_Float_Matrix(x_size, y_size);
   out_img = Alloc_Short_Matrix(x_size, y_size);
   char_img = Alloc_Char_Matrix(x_size, y_size);
   File2_Char_Matrix(in_f, char_img, x_size, y_size);
   for(i=0;i<x_size;i++)
   {	   
     for(j=0;j<y_size;j++)
	 {
		 in_img[i][j]=char_img[i][j];
		 printf("%4d ",char_img[i][j]);
	 }
	 printf("\n");
   }

   printf("\n");

   f97_2D(in_img, x_size, y_size, levels, 0);
    for(i=0;i<x_size;i++)
	{
     for(j=0;j<y_size;j++)
	 {
       out_img[i][j]=(short)in_img[i][j];
	   printf("%4d ",out_img[i][j]);
	 }
	 printf("\n");
	}
	printf("\n");
   Short_Matrix_2File(out1_f, out_img, x_size, y_size);


   f97_2D(in_img, x_size, y_size, levels, 1);
   for(i=0;i<x_size;i++)
     for(j=0;j<y_size;j++)
      {
       if(in_img[i][j]>255)
          char_img[i][j]=255;
          else if(in_img[i][j]<0)
            char_img[i][j]=0;
            else
             char_img[i][j]=in_img[i][j];
      }
	 
   for(i=0;i<x_size;i++)
	{
     for(j=0;j<y_size;j++)
	 {       
	   printf("%4d ",char_img[i][j]);
	 }
	 printf("\n");
	}

   Char_Matrix_2File(out2_f, char_img, x_size, y_size);

  printf("psnr=%f \n",PSNR(in_f,out2_f,x_size,y_size));
  */
   readfile("F:\\MATLAB71\\work\\s64.wav",&sample_data,&datalen);//framelen 数据总长度（字节）；

   printf("sample_data=%4X\n",sample_data);
   printf("datalen=%d\n",datalen);
  for(i=0;i<datalen;i++)
  {    
	   data_f[i]=(float)sample_data[i];
  }
   f97_1D(data_f, datalen, 2, 0);
   f97_1D(data_f, datalen, 2, 1);

  for(i=0;i<8;i++)
  {    
	   printf("%4f ",x[i]);
  }
  printf("\n");
  f97_1D(x, 8, 1, 0);
  for(i=0;i<8;i++)
  {    
	   printf("%4f ",x[i]);
  }
  printf("\n");

  f97_1D(x, 8, 1, 1);
  for(i=0;i<8;i++)
  {    
	   printf("%4f ",x[i]);
  }
  printf("\n");
  
  //getchar()
  while(1)
  {}

}