维纳滤波器设计

1.设计要求

Sequence s(n) of N=2000 points is generated by AR(1) model: s(n)=as(n-1)+w(n), in which a=0.8, w(n) is white noise sequence, the mean and variance of w(n) ismw0，2w0.36.

The measurement model is x(n) =s(n) +v(n), in which white noise sequence v (n) and

21. w (n) is not related, the mean and variance of v(n) is mv0,m

Requirements:

(1)Design IIR causal Wiener filter , calculate the filtered sequence and mean square error;

(2)Design FIR Wiener filter , calculate the filtered sequence and mean square error;

(3)Display raw data , noise data and filtered data on the same graph , compare the mean square error between the two cases and draw a conclusion.

2.设计原理

2.1维纳滤波原理概述

维纳（Wiener）是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。一个线性系统，如果它的单位样本响应为h(n)，当输入一个随机信号x(n)，且

x(n)s(n)v(n) （1）其中x(n)表示信号，v(n))表示噪声，则输出y(n)为

y(n) h(m)x(nm) （2）

我们希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因此称y(n)为s(n)的估计值，用s(n)表示，即

y(n)s(n) （3） ^^

则维纳滤波器的输入—输出关系可用下面图1表示。

图1

实际上，式（2）所示的卷积形式可以理解为从当前和过去的观察值x(n)，x(n1)，

x(n2)„x(nm)，„来估计信号的当前值s(n)。因此，用h(n)进行过滤问题^

实际上是一种统计估计问题。

一般地，从当前的和过去的观察值x(n)，x(n1)，x(n2)„估计当前的信号值y(n)s(n)成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值y(n)s(nN)(N0)称为外推或预测；从过去的观察值，估计过去的信号值y(n)s(nN)(N1)称为平滑或内插。因此维纳滤波器又常常被称为最佳线性^^^

过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。

如果我们分别以s(n)与s(n)表示信号的真实值与估计值，而用e(n)表示他们之间的误差，即

e(n)s(n)s(n) （4）显然e(n)可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小：

(n)E[ e(n)]min （5）2^^

采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。

2.2维纳-霍夫方程的求解

为了按（5）式所示的最小均方误差准则来确定维纳滤波器的冲激响应h(n)，令(n)对h(j)的导数等于零，即可得

Rxs(m)h(i)R

ixx (mi),m （6）

式中，Rxs(m)是s(n)与x(n)的互相关函数，Rxx(m)是x(n)的自相关函数，分别定义为

RxsE[x(n)s(nm)]

RxxE[x(n)x(nm)]

式（6）称为维纳滤波器的标准方程或维纳-霍夫（Wiener-Hopf）方程。如果已知Rxs(m)和Rxx(m)，那么解此方程即可求的维纳滤波器的冲激响应。

式（6）所示标准方程右端的求和范围即i的取值范围没有具体标明，实际上

有三种情况：

（1）有限冲激响应（FIR）维纳滤波器，i从0到N1取得有限个整数值；

（2）非因果无限冲激响应（非因果IIR）维纳滤波器，i从到取所

有整数值；

（3）因果无限冲激响应（因果IIR）维纳滤波器，i从0到取正整数值。上述三种情况下标准方程的解法不同，本文将描述因果IIR维纳滤波器和FIR维纳滤波器的求解。

3维纳滤波器的设计与实现

3.1因果IIR维纳滤波器的求解

解维纳-霍夫方程得，因果IIR维纳滤波器的传输函数为：

H1

cz2Ssxz

1

BzBz

现要设计一个因果IIR维纳滤波器，对x(n)进行处理，以得最佳估计。

为分析简单起见，作如下假设：

（1）wn是方差等于Q的白噪声，其自相关函数和功率谱分别为

EwnwiQ1,ni

ni ni0,ni

（2）vn是方差等于R的白噪声，其自相关函数和功率谱分别为

EvnviRni，

且 SwwzQ

SvvzR

（3）vn与s(n)不相关，与wn也不相关，即

Evnsi0

E

vnwi0

我们需要求出如下几个功率谱：

SsszQ1az11az

ScQ

sxz1az11az 7）（

c2QSxxzcSsszSvvzR 11az1az2

对出入信号功率谱进行因式分解，得

1fz1

SxxzBzBz,Bz,1az121f1 （8）

继续推导得到

21f2cQ1a2R

（9） faR1f2c2Q1a2R

联立两个方程，得

a2RPQP,P0 （10） 2RcP

求得因果IIR滤波器系统函数

Hcz1Gz 2Bz1cQ11 （11） G1111fa1fz1az1fz21az1

其中，维纳增益 G

fcP （12） 2RcPaRa1cG （13） Rc2P

3.2 FIR维纳滤波器的求解

设滤波器冲激响应序列的长度为N，冲激响应矢量为

h[h(0)h(1)....h(N1)] （14） T

滤波器输入数据矢量为

x(n)[x(n)x(n1)...x(nN1)] （15） T

则滤波器的输出为

TT y(n)s(n)x(n)hhx(n) （16） ^

这样，式（6）所示的维纳-霍夫方程可写成

PhR或PRh （17） TT

其中

PE[x(n)s(n)] （18）是s(n)与x(n)的互相关函数，它是一个N维列矢量；R是x(n)的自相关函数，是N阶方阵

RE[x(n)x(n)] （19） T

利用求逆矩阵的方法直接求解式（10），得

hoptRP （20） 1

这里opt表示“最佳”，这就是FIR维纳滤波器的冲激响应。

3.3 维纳滤波器的matlab实现

由设计要求，可得维纳滤波器的数据模型和测量模型参数：

a = 0.8 , c = 1 , Q = 0.36 , R = 1 ;

由3.1节，我们可得IIR维纳滤波器的设计步骤如下：

（1）根据一直参数a , c , R , Q求解式（10），得到正解P = 0.6 ;

（2）由式（12）计算维纳增益G= 0.375；

（3）由式（13）计算滤波器系数f= 0.5；

（4）将G和f值带入式（11），得到IIR维纳滤波器的传输函数 Hcz

于是单位脉冲响应为h(n)0.375*(0.5)nu(n)。

而由3.2节的计算推导，我们可以设计FIR 维纳滤波器的程序流程图为：

0.375，10.5z1

图1 FIR维纳滤波器算法流程图

RE[x(n)x(n)] （19） T

利用求逆矩阵的方法直接求解式（10），得

hoptRP （20） 1

这里opt表示“最佳”，这就是FIR维纳滤波器的冲激响应。

3.3 维纳滤波器的matlab实现

由设计要求，可得维纳滤波器的数据模型和测量模型参数：

a = 0.8 , c = 1 , Q = 0.36 , R = 1 ;

由3.1节，我们可得IIR维纳滤波器的设计步骤如下：

（1）根据一直参数a , c , R , Q求解式（10），得到正解P = 0.6 ;

（2）由式（12）计算维纳增益G= 0.375；

（3）由式（13）计算滤波器系数f= 0.5；

（4）将G和f值带入式（11），得到IIR维纳滤波器的传输函数 Hcz

于是单位脉冲响应为h(n)0.375*(0.5)nu(n)。

而由3.2节的计算推导，我们可以设计FIR 维纳滤波器的程序流程图为：

0.375，10.5z1

图1 FIR维纳滤波器算法流程图

由此，具体Matlab代码如下：

clc

clear all

L = 2000; % signal length

N = 20; % length of the FIR filter

a = 0.8;

% white noise with mean of 0 and var of 0.36

w = sqrt( 12 * 0.36) * ( rand(1,L) - 0.5 );

% true signal: s(n)=a*s(n-1)+w(n)

s = zeros(1,L); s(1) = w(1);

for ii = 2:L

s(ii) = a * s(ii-1) + w(ii);

end

% white noise with mean of 0 and var of 1

v = sqrt( 12 ) * ( rand(1,L) - 0.5 );

% received signal: x = s + v

x = s + v;

% r_xx is the autocorrelation of x

r_xx = xcorr( x );

% R_xx is the N-dimentional autocorrelation matrix of x R_xx = zeros( N );

for ii = 1 : N

R_xx( : , ii ) = r_xx( L+1-ii : L+1-ii+N-1 )';

end

% r_xs is the cross-correlation of x and s

r_xs = xcorr( x , s );

r_xs = r_xs( L : L+N-1 )';

% according to R_xx * h_FIR = R_xs ,we can calculate the h_FIR h_FIR = inv(R_xx) * r_xs;

h_FIR = h_FIR';

% y_FIR is the reslut signal filtered by h_FIR

y_FIR = cconv( h_FIR , x, L );

% y_IIR is the result signal filtered by IIR causal Wiener filter h_IIR

n = 0:L-1;

h_IIR = 0.375 * ( 0.5 .^ n );

y_IIR = cconv( h_IIR , x, L );

%% plot the true signal s(n) , the received signal x(n), the result signal y_FIR(n)

%%filtered by h_FIR,the result signal y_IIR(n) filtered by h_IIR of the last 100 points

t = L - 99 : L;

figure(1);

subplot(2,2,1);

plot( t , x(t), 'k' );

title('The raw data s(n) ');