lp.c

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

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

#define PI 3.14159265
#define FRAMES 1
#define FRAME_LEN 80
#define ORDER 20

int ReadFrame(FILE * data, signed short *);
int Windowing(signed short *);
int EstimateAutoCorrelation(signed short *, signed int *);
double * Levinson_Durbin(signed int *, int);
int Prediction(signed short *, double *, signed short *);
int Error(signed short *, signed short *, signed short *, signed short *);
int PrintData(signed short *);
int PrintR(signed int *);
int PrintA(double *);

int main()
{
	signed short frame_data[FRAME_LEN];
	signed int r[FRAME_LEN];
	signed short prediction_data[FRAME_LEN];
	signed short error[FRAME_LEN];
	double * a;
	unsigned short se;
	int i;
	int frame = 0;
	FILE * data = fopen("u.dat", "rb");
	FILE * SOURCE = fopen("source.txt", "w");
	FILE * PREDICT = fopen("predict.txt", "w");
	
	if(data == NULL){
		printf("Error: Unable to open data file!\n");
		return 1;
	}
	if(SOURCE == NULL){
		printf("Error: Unable to open source file!\n");
		return 1;
	}
	if(PREDICT == NULL){
		printf("Error: Unable to open predict file!\n");
		return 1;
	}

	fprintf(SOURCE, "s=[");
	fprintf(PREDICT, "p=[");
	while(!ReadFrame(data, frame_data)){
	  for(i=0;i<FRAME_LEN;i++)
	  {
		  fprintf(SOURCE, "%d ", frame_data[i]);
	  }
	  EstimateAutoCorrelation(frame_data, r);
	  a = Levinson_Durbin(r, ORDER);
      Prediction(frame_data, a, prediction_data);
	  Error(frame_data, prediction_data, error, &se);
	  for(i=0;i<FRAME_LEN;i++)
	  {
		  fprintf(PREDICT, "%d ", prediction_data[i]);
	  }
	  printf("Frame %d MSE=%u\n", frame+1, se);
	  if(++frame>=FRAMES)
	  {
		  break;
	  }
	}	
	fprintf(SOURCE, "];");
	fprintf(PREDICT, "];");

	printf("All data processed!\n");

	fclose(data);
	fclose(SOURCE);
	fclose(PREDICT);

	return 0;
}

int ReadFrame(FILE * data, signed short * FramePtr)
{
	if(fread(FramePtr, sizeof(signed short), FRAME_LEN, data) < FRAME_LEN){
		return 1;
	}else{
		return 0;
	}
}

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

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

  return 0;
}

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

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

	return 0;
}

double * Levinson_Durbin(signed int * RPtr, int M)
{
	int i,j;
	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;
        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);			
	}

	//printf("%f\n", p);
	return a;
}

int Prediction(signed short * FramePtr, double * a, signed short * PredictionPtr)
{
	int i,j;
	static buffer[ORDER] = {0};
	double tmp, sum;

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

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

	tmp = (double)(FRAME_LEN);
	sum = 0;
	for(i=0;i<FRAME_LEN;i++)
	{
		err[i] = source[i] - predict[i];
		//if(i>=ORDER)
		//{
		  sum += ((double)err[i] * (double)err[i])/tmp;
		//}
	}	
	//printf("%f\n", sum);
	*se = (signed short)sqrt(sum);

	return 0;
}

int PrintData(signed short * data)
{
	int i;

	for(i=0;i<FRAME_LEN;i++)
	{
		printf("%04X ", data[i]);
		if(i%8==7)
		{
			printf("\n");
		}
	}

	return 0;
}

int PrintR(signed int * r)
{
	int i;
	for(i=0;i<FRAME_LEN;i++)
	{
		printf("%08X ", r[i]);
		if(i%8==7)
		{
			printf("\n");
		}
	}

	return 0;
}

int PrintA(double * a)
{
	int i;

	for(i=0;i<=ORDER;i++)
	{
		printf("%8f\t", a[i]);
		if(i%5==4)
		{
			printf("\n");
		}
	}

	return 0;
}