burg.c

来自「DSP信号处理源码,包括数字信号处理课程中的基本源程序。」· C语言 代码 · 共 91 行

/**************************************************************
Burg 算法
***************************************************************/
#include <stdio.h>
#include <math.h>
#include <graphics.h>
#include <stdlib.h>

#define P 5                                       /*阶数*/
#define LENGTH 15                                 /*样点长度*/
#define F     20                                  /*取样频率*/
#define F1    2                                   /*信号频率*/
#define F2    8                                   /*信号频率*/

/**************************************************************
*  int    n         :输入数据的个数                           *
*  int    buffer[]  :输入数据缓冲区，用于存放被分析的整数数据 *
*  floar  a[p+1]    :用于存放计算出的系数                     *
***************************************************************/
float       burg(int n, int *buffur, float *a);
main()
{
int             i,buffer[LENGTH];
float           a[P+1];
float           ew,f1=F1,f2=F2,f=F;
int             p=P,length=LENGTH;
   for(i=0;i<LENGTH;i++)
        {
	  buffer[i]=(int)1000*(sin(F1*2*3.1415926*i/F)+
	  sin(F2*2*3.1415927*i/F));
	}

  ew=burg(LENGTH,buffer,a); getch();
}


float burg(int n,int *buffur,float *a)                   /*Burg算法*/
{  int      i,m;
  float    ew=0;
  float    pm,pm1,pm2;
  float    ate[P+1], k[P+1], ef[LENGTH],
  eb[LENGTH],efte[LENGTH],ebte[LENGTH];
  a[0]=1;                                            /*初始化数据*/
  for(i=0;i<n;i++)
      {
	  ef[i]=eb[i]=buffur[i];
	  ew=ew+(float)((float)(1/(float)(n))*(float)(buffur[i])*(float)
	  (buffur[i]));
      }
  for(m=1;m<=P;m++)                                  /*计算模型系数*/
      {
	   pm1=pm2=0;

	   for(i=m;i<n;i++)                          /*计算反射系数*/
	   {
	       pm1+=ef[i]*eb[i-1];
	       pm2+=ef[i]*ef[i]+eb[i]*eb[i];
	   }

	   pm=-2*pm1/pm2;

	   for(i=0;i<=m;i++)                          /*计算滤波器系数*/
	   {
	       ate[i]=a[i];
	   }

	   a[m]=pm;

	   for(i=1;i<m;i++)
	   {
	       a[i]=ate[i]+pm*ate[m-i];
	   }

	   ew=(1-pm*pm)*ew;                           /*计算预测误差功率*/

	   for(i=0;i<n;i++)                           /*计算前后向预测误差*/
	   {
	       efte[i]=ef[i];
	       ebte[i]=eb[i];
	    }

	   for(i=0;i<n;i++)
	   {
	       ef[i]=efte[i]+pm*ebte[i-1];
	       eb[i]=ebte[i-1]+pm*efte[i];
	    }
      }
      return(ew);
}