FIR维纳滤波的设计

FIR维纳滤波器的设计在信号处理的许多实际应用中，人们往往不能直接获得所需要的有用信号，需要从噪声中提取信号。比如，在信号传输过程中，由于存在信道噪声等干扰，在接收端观测到的信号必然与原始信号不同。为了从观测数据中尽可能精确地重现原始信号，而最大成都地抑制噪声，需要设计一种滤波器，其输出尽可能逼近原始信号，成为原始信号的最佳估计。这种滤波器成为最佳滤波器。维纳(Wiener)滤波器就是用来解决这样一类问题的一种滤波器。本文将应用MATLAB并结合实例介绍FIR维纳滤波器的设计方法。一、维纳滤波的原理维纳滤波的本质是一种最佳估计问题，采用的是最小均方误差准则。一个线性系统，其单位样本响应为h(n)，当输入一个随机信号其中s(n)表示信号，表示噪声，则输出y(n)为 (1)系统是通过y(n)来估计s(n)，因此将其称为s(n)的估计值，用表示，即 (2)图1维纳滤波器基本框图图1所示为维纳滤波器的基本框图。式(1)为一卷积，可以理解为从当前和过去的观察值x(n)，x(n-1)，x(n-2)…x(n-m)，…来估计信号的当前值。维纳滤波器一般有三种用途。用当前的和过去的观察值x(n)，x(n-1)，x(n-2)，…来估计当前的信号值称为滤波；用过去的观察值来估计当前的或将来的信号值称为预测；用全部数据来估计过去的信号值称为平滑。维拉滤波采用的是最小均方误差准则，以与分别表示信号的真值与估计值，而用e(n)表示它们之间的误差，即 (3)那么最小均方误差表示为 = (4)为了得到使式(4)最小时的维纳滤波器冲击响应，令对的导数等于零，即(5)由式(5)可得(6)式(6)称为正交方程，表明任何时刻的估计误差都与用于估计的所有数据正交。将式(1)，式(3)代入式(6)可得(7)式(7)称为维纳-霍夫（Wiener-Hopf）方程，其中是与的互相过函数，是的自相关函数，其表达式如下(8)(9)维纳滤波器有三种情况，一是FIR维纳滤波器；二是非因果IIR维纳滤波器；三是因果的IIR维纳滤波器。本文仅讨论FIR维纳滤波器的求解。二、FIR维纳滤波器求解设维纳滤波器的冲击响应序列长度为N，则冲击响应矢量为h=(10)滤波器输入数据矢量为=(11)着滤波器输出为(12)这样，维纳-霍夫方程可表示成P=Rh(13)其中P=，R=。求解式(13)可得h=三、实例设计 结合维拉滤波器的原理，利用MATLAB设计滤波器。 将随机信号X(n)看成是由典型白噪声序列源W(n)激励一个线性系统产生，用一个差分方程来描述。进行Z变换得到，那么均值为１的高斯白噪声序列W（n）可以用randn函数产生，再利用函数X=filter(B,A,W)产生随机信号x(n)。 这里将滤波器的阶数设为101，根据维纳－霍夫方程：，其中是观测信号的自相关函数，是观测信号和期望信号的互相关函数。定义维纳滤波的模型，最后带入filter。具体MATLAB程序代码如下：maxlag=100;N=100; %采样次数为100x=zeros(N,1);y=zeros(N,1);var=1; %%%%%%%%%%%%%%%%%%列出状态方程x(1)=randn(1,1);%令x(-1)=x(-2)=x(-3)=x(-4)=0x(2)=randn(1,1)+1.352*x(1);x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);forn=5:Nx(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);end;v=randn(N,1);y=x+v; %z_x为观测样本值=真值+噪声%%%%%%%%%%%%%%%%%%%%%%%滤波x=x';y=y';xk_s(1)=y(1); %赋初值xk_s(2)=y(2);xk_s(3)=y(3);xk_s(4)=y(4);xk=[y(1);y(2);y(3);y(4)];%%%%%%%%%%%%%%%%%%计算观测函数%维纳滤波器的生成[rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵，滤波器的阶数为101阶rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数rx2=rx2(101:end);h=inv(rx1)*rx2'; %维纳-霍夫方程xk_s=filter(h,1,y);%加噪信号通过滤波器后的输出%%%%%%%%%%%%%%%%%%%%%计算误差e_x=0;eq_x=0;e_x1=N:1;%计算滤波的均值,计算滤波误差的均值fori=1:Ne_x(i)=x(i)-xk_s(i);%误差=真实值-滤波估计值end%%%%%%%%%%%%%%%%%%%%%%%作图t=1:N;figure(1);plot(t,x,'r',t,y,'g',t,xk_s,'b');legend('