lp2.c

自适应滤波中线性预测算法Visual C++ 6.0代码

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

#define PI 3.14159265

int ReadFrame(FILE * data, signed short *);     //read a frmae of data from the input file
int Windowing(signed short *, signed short *);  //windowing the frame
int EstimateAutoCorrelation(signed short *, signed int *);  //estimate autocorrelation
double * Levinson_Durbin(signed int *, int);                //Levninson-Durbin Algorithm
int Prediction(signed short *, double *, signed short *);   //Prediction
int Error(signed short *, signed short *, signed short *, double *);   //calculate prediction error

int ORDER;                  //prediction order
int WINDOW_LEN;             //window length
int FRAMES;                 //frames to be processed
int FRAME_LEN;              //frame length

int main()
{
	signed short * frame_data;           //frame data
	signed short * window_data;          //frame data after windowing
    signed int * r;                      //autocorrelatin function
	signed short * prediction_data;      //prediction data
	signed short * error;                //prediction error
	double * a;                          //filter coefficients
	double se = 0;                       //square error
	int exit = 0;
	int i;
	int frame = 0;
	int start, stop, step;
	int window;                          //1:Rectangular window 2:Hanning window
	FILE * data;                         //input data file
	FILE * SOURCE = fopen("source.txt", "w");
	FILE * COEFFICIENT = fopen("coefficient.txt", "w");
	FILE * PREDICT = fopen("predict.txt", "w");
	FILE * DIFF = fopen("diff.txt", "w");

	if(SOURCE == NULL){
		printf("Error: Unable to open source file!\n");
		return 1;
	}
	if(COEFFICIENT == NULL){
		printf("Error: Unable to open coeffiecient file!\n");
		return 1;
	}
	if(PREDICT == NULL){
		printf("Error: Unable to open predict file!\n");
		return 1;
	}
	if(DIFF == NULL){
		printf("Error: Unable to open diff file!\n");
		return 1;
	}	

	fprintf(SOURCE, "s=[");
    fprintf(COEFFICIENT, "a=[");
	fprintf(PREDICT, "p=[");
	fprintf(DIFF, "d=[");

	printf("Input Frame Length:\n");
	printf("?:");
	scanf("%d", &FRAME_LEN);
	if(FRAME_LEN<=0)
	{
		printf("Error: Wrong Input!\n");
		return 1;
	}

	printf("Input frames to be processed:\n");
	printf("?:");
	scanf("%d", &FRAMES);
	if(FRAMES<=0)
	{
		printf("Error: Wrong Input!\n");
		return 1;
	}	

	printf("Select Window:\n1 Rectangle\n2 Hanning\n?:");
	scanf("%d", &window);
	if(window != 1 && window != 2)
	{
		printf("Wrong Input!\n");
		return 1;
	}
	printf("Input start order, stop order and step:\n");
	scanf("%d%d%d", &start, &stop, &step);
	if(start <= 0 || stop >= FRAME_LEN || stop < start)
	{
		printf("Wrong Input!\n");
		return 1;
	}

	prediction_data = malloc(FRAME_LEN*sizeof(signed short));
    error = malloc(FRAME_LEN*sizeof(signed short));

	for(ORDER=start;ORDER<=stop;ORDER+=step)
	{	  
	  WINDOW_LEN = FRAME_LEN + ORDER;
	  frame_data = malloc(WINDOW_LEN*sizeof(unsigned short));
      window_data = malloc(WINDOW_LEN*sizeof(unsigned short));
	  r = malloc((ORDER+1)*sizeof(signed int));
      for(i=0;i<ORDER;i++)
	  {
		frame_data[i] = 0;
	  }
	  se = 0;
	  frame = 0;
	  data = fopen("u.dat", "rb");
	  if(data == NULL){
	    printf("Error: Unable to open data file!\n");
		return 1;
	  }
	  while(!ReadFrame(data, &frame_data[ORDER])){        //read data from input file
	    for(i=ORDER;i<WINDOW_LEN;i++)
		{
		  fprintf(SOURCE, "%d ", frame_data[i]);
		}
		if(window == 1)
		{
          EstimateAutoCorrelation(frame_data, r);         
		}else{
	      Windowing(frame_data, window_data);             //windowing
	      EstimateAutoCorrelation(window_data, r);        //estimate autocorrelation
		}
	    a = Levinson_Durbin(r, ORDER);                    //L-D Algorithm
        Prediction(frame_data, a, prediction_data);       //prediction
		for(i=0;i<=ORDER;i++)
		{
			fprintf(COEFFICIENT, "%f ", a[i]);
		}
	    Error(&frame_data[ORDER], prediction_data, error, &se);    //calculate prediction error
	    for(i=0;i<FRAME_LEN;i++)
		{
		  fprintf(PREDICT, "%d ", prediction_data[i]);
		}
	    for(i=0;i<FRAME_LEN;i++)
		{
		  fprintf(DIFF, "%d ", error[i]);
		}
	    for(i=0;i<ORDER;i++)
		{
		  frame_data[i] = frame_data[WINDOW_LEN-ORDER+i];
		}	 
	    if(++frame>=FRAMES)
		{
		  break;
		}
	  }	

	  if(ORDER+step<=stop){
	    fprintf(SOURCE, ";");
		fprintf(COEFFICIENT, ";");
	    fprintf(PREDICT, ";");
	    fprintf(DIFF, ";");
	  }      
	  printf("ORDER=%-5d MSE=%.0f\n", ORDER, se);
	  free(frame_data);
	  free(window_data);
	  free(r);
	  fclose(data);
	}

	fprintf(SOURCE, "];\n");
	fprintf(COEFFICIENT, "];\n");
	fprintf(PREDICT, "];\n");
	fprintf(DIFF, "];\n");
	fclose(SOURCE);
	fclose(PREDICT);
	fclose(DIFF);

	printf("Process Complete!\nWould you like to exit now(1 to exit)?\n");
	scanf("%d", &exit);
	while(!exit)
	{
		printf("?:");
        scanf("%d", &exit);
	}

	return 0;
}

/*Read a frame of data from the input file*/
int ReadFrame(FILE * data, signed short * FramePtr)
{
	if(fread(FramePtr, sizeof(signed short), FRAME_LEN, data) < (unsigned)FRAME_LEN){
		return 1;
	}else{
		return 0;
	}
}

/*Hanning Windowing*/
int Windowing(signed short * FramePtr, signed short * WindowPtr)
{
  int i;
  double tmp1, tmp2;

  tmp1 = 2*PI/(WINDOW_LEN - 1);
  for(i=0;i<WINDOW_LEN;i++)
  {
	  tmp2 = 0.5 - 0.5*cos(i*tmp1);
	  tmp2 = tmp2 * FramePtr[i];
      WindowPtr[i] = (signed short)tmp2;
  }

  return 0;
}

/*Estimate autocorrelation function*/
int EstimateAutoCorrelation(signed short * FramePtr, signed int * RPtr)
{
	int i, j;	
	double tmp, sum;

	for(i=0;i<ORDER+1;i++)
	{
		sum = 0;
		for(j=0;j<WINDOW_LEN-i;j++)
		{
			tmp = (double)FramePtr[j] * (double)FramePtr[j+i];				
			tmp = tmp / WINDOW_LEN;
			sum += tmp;
		}		
		RPtr[i] = (signed int)sum;
	}

	return 0;
}

/*Levinson-Durbin Algorithm*/
double * Levinson_Durbin(signed int * RPtr, int M)
{
	int i,j;
	int warning = 0;
	double delta = (double)RPtr[1];
	double p = (double)RPtr[0];
	double k;
	double * a = malloc(M*sizeof(double));
	double * tmp = malloc(M*sizeof(double));
	a[0] = 1.0;
	tmp[0] = 0;

	for(i=1;i<=M;i++)
	{
		delta = 0;
		for(j=0;j<i;j++)
		{
			delta += (double)RPtr[i-j] * a[j];
		}
		k = -delta/p;
		if(k<-1.0 || k>1.0)
		{
			warning = 1;			
		}
        a[i] = 0;		
		for(j=1;j<=i;j++)
		{
			tmp[j] = a[i-j];
		}
		for(j=0;j<=i;j++)
		{
			a[j] = a[j] + k*tmp[j];
		}
		p = p*(1-k*k);			
	}

	if(warning)
	{
	  printf("WARNING: |k| > 1 when ORDER=%d The error-filter will be unstable!\n", ORDER);
	}

	return a;
}

/*Linear Prediction*/
int Prediction(signed short * FramePtr, double * a, signed short * PredictionPtr)
{
	int i,j;
	double tmp, sum;

	for(i=0;i<FRAME_LEN;i++)
	{
       sum = 0;
	   for(j=1;j<=ORDER;j++)
	   {
		   tmp = -a[j] * FramePtr[ORDER+i-j];
		   sum += tmp;
	   }
       PredictionPtr[i] = (signed short)sum;
	}

	return 0;
}

/*Calculate Prediction Error*/
int Error(signed short * source, signed short * predict, signed short * err, double * se)
{
	int i;
	double tmp, sum;

	tmp = (double)(FRAME_LEN*FRAMES);
	sum = 0;
	for(i=0;i<FRAME_LEN;i++)
	{
		err[i] = source[i] - predict[i];
		sum += ((double)err[i] * (double)err[i])/tmp;		
	}	
	*se += sum;

	return 0;
}