wienerexam.c

自适应维纳滤波器函数及调用方法

//---------------------------------------------------------------------------

#pragma hdrstop

//---------------------------------------------------------------------------

#pragma argsused
#include "stdio.h"
#include "math.h"
void wiener(double rxx[],double rdx[],int p,double h[],double *e,
      double x[],double y[],int n);
int levin(double t[],double b[],int n,double x[]);
void main()
{
  int i,p,n;
  double e,rxx[3],rdx[3],h[3],x[10],y[10];
  n=1;
  p=2;
  for(i=0;i<=p;i++)
  {rdx[i]=4.0*pow(0.5,i);}
  rxx[0]=rdx[0]+2.0;
  for(i=1;i<=p;i++)
  {rxx[i]=rdx[i];}
  wiener(rxx,rdx,p,h,&e,x,y,n);
  printf("\n The Coefficients of Wiener Filter\n");
  for(i=0;i<=p;i++)
    {printf("h(%d)=%10.7lf\n",i,h[i]);}
  printf("The Minimum MSE Error=%lf\n",e);
  while(getchar()==' ');
}

/* rxx[p+1],信号x（n）的自相关函数rxx（i）
   rdx[p+1]，信号d（n）与x（n）的互相关函数rdx（i）
   p 整型变量，维纳滤波器的阶数
   h[p+1]，微纳滤波器的系数。
   e 双精度实型变量。维纳滤波器的最小均方误差
   x[n] 输入信号
   y[n] 维纳滤波后的输出信号
   n 输入信号长度
*/
void wiener(double rxx[],double rdx[],int p,double h[],double *e,
      double x[],double y[],int n)
{
  int i,k;
  double sum;
  levin(rxx,rdx,p+1,h);
  sum=0.0;
  for(i=0;i<=p;i++)
  {
     sum+=rdx[i]*h[i];
  }
  *e=rdx[0]-sum;
  for(k=0;k<n;k++)
  {
     y[k]=0.0;
     for(i=0;i<=p;i++)
     {
        if((k-i)>=0)y[k]+=h[i]*x[k-i];
     }
  }
}

/* t[n] 存放对称托步利兹矩阵的元素t0，t1，t2，...，tn-2，tn-1
   b[n] 方程组右端的常数向量
   n 方程组的阶数
   x[n] 方程组的解
*/
int levin(double t[],double b[],int n,double x[])
{
   int i,j,k;
   double a,beta,q,c,h,*y,*s;
   s=malloc(n*sizeof(double));
   y=malloc(n*sizeof(double));
   a=t[0];
   if((fabs(a)+1.0)==1.0)
   {
      free(s);
      free(y);
      printf("ill_conditioned\n");
      return(-1);
   }
   y[0]=1.0;
   x[0]=b[0]/a;
   for(k=1;k<n;k++)
   {
     beta=0.0;
     q=0.0;
     for(j=0;j<k;j++)
     {
        beta+=y[j]*t[j+1];
        q+=x[j]*t[k-j];
     }
   if((fabs(a)+1.0)==1.0)
   {
      free(s);
      free(y);
      printf("ill_conditioned\n");
      return(-1);
   }
   c=-beta/a;
   s[0]=c*y[k-1];
   y[k]=y[k-1];
   if(k!=1)
   {
      for(i=1;i<k;i++)
      {s[i]=y[i-1]+c*y[k-i-1];}
   }
   s[k]=y[k-1];
   a+=c*beta;
   if((fabs(a)+1.0)==1.0)
   {
      free(s);
      free(y);
      printf("ill_conditioned\n");
      return(-1);
   }
   h=(b[k]-q)/a;
   for(i=0;i<k;i++)
   {
     x[i]+=h*s[i];
     y[i]=s[i];
   }
   x[k]=h*y[k];
 }
 free(s);
 free(y);
 return(1);

}