实验一维纳滤波器的计算机实现.docx

实验一维纳滤波器的计算机实现.docx

《实验一维纳滤波器的计算机实现.docx》由会员分享，可在线阅读，更多相关《实验一维纳滤波器的计算机实现.docx（16页珍藏版）》请在冰豆网上搜索。

实验一维纳滤波器的计算机实现

实验一维纳滤波器的计算机实现

一、设计目的

1．利用计算机编程实现加性干扰信号的维纳滤波。

2．将计算机模拟实验结果与理论分析结果相比较，分析影响维纳滤波效果的各种因素，从而加深对维纳滤波的理解。

二、设计原理与方法

维纳滤波是一种从噪声背景中提取信号的最佳线性方法，假定一个随机信号x（n）具有以下形式：

（1-1）

其中，s（n）为有用信号，v（n）为噪声干扰，将其输入一个单位脉冲响应为h（n）的线性系统，其输出为

（1-2）

我们希望x（n）通过这个系统后得到的y（n）尽可能接近于s（n），因此，称y（n）为信号s（n）的估值。

按照最小均方误差准则，h（n）应满足下面的正则方程：

（1-3）

这就是著名的维纳－霍夫方程，其中

（1-4）

是x（n）与s（n）的互相关函数，定义为

（1-5）

这里，E[·]表示求数学期望，*表示取共轭。

在要求h（n）满足因果性的条件下，求解维纳-霍夫方程是一个典型的难题。

虽然目前有几种求解h（n）的解析方法，但它们在计算机上实现起来非常困难。

因此本实验中利用近似方法，即最佳FIR维纳滤波方法，在计算机上实现随机信号的维纳滤波。

设h（n）为一因果序列，其长度为N，则

（1-6）

同样利用最小均方误差准则，h（n）满足下面正则方程：

（1-7）

其中

（1-8）

（1-9）

（1-10）

这里T表示转置运算。

称

为信号x（n）的N阶自相关矩阵，

为x（n）与s（n）的互相关函数向量。

当

为满秩矩阵时，由公式（1-7）可得

（1-11）

由此可见，利用有限长的h（n）实现维纳滤波器，只要已知

和

，就可以按上式解得满足因果性的h。

虽然它不同于真正的维纳滤波器，但是只要N选择的足够大，它就可以很好地逼近真正的维纳滤波器，这一点我们可以在下面实验中得到证实。

在本实验中，s（n）由下式来确定：

（1-12）

其中

是零均值方差为

的均匀分布白噪声，v（n）是与s（n）互不相关的均匀分布白噪声，其均值为零，方差

。

根据理论推导，此时维纳最佳滤波器为

（1-13）

单位脉冲响应为

（1-14）

由此可以实现对信号x（n）的最佳过滤，即

（1-15）

其中

为s（n）的最佳估值。

同时可以推出，经过理想维纳滤波后，均方误差应为

（1-16）

在实验中，我们利用下面公式来统计均方误差：

（1-17）

其中L为维纳滤波数据长度。

实际中，一般很难确切地知道

和

，通常是利用有限个x（n）和s（n）的样本来估计它们

（1-18）

（1-19）

为了在检验实际中某次产生序列的自相关特性与理论值的近似程度，我们可以采用下式进行度量：

（1-20a）

（1-20b）

该式表示了自相关函数的理论值与某次实现的实际值的相对平方误差。

实验中为了得到与自相关特性理论值相符的观测序列，往往需要多次产生序列，直到两者的相对平方误差

足够小。

本实验中，我们取K=50，并认为

且

的序列才是满足要求的。

三、设计内容与步骤

1．仔细阅读维纳滤波原理，根据图1.1给出的框图编制维纳滤波程序。

2．运行维纳滤波程序，选择L=500,N=10，观察并记录实验结果，分析比较下列三个问题：

①与s（n）比较，信号x（n）在维纳滤波前后有何差别？

滤波效果如何？

②估计出

的和理想的h（n）比较，近似程度如何？

③理想的维纳滤波和FIR维纳滤波效果有何差异？

图1.1

3．固定L=500,分别取N=3、20,根据实验结果，观察N的大小对

的估计和滤波效果的影响。

记录实验结果。

4．固定N=10,改变L=200,1000,根据实验结果，观察并记录L的大小对

的精度和滤波效果的影响。

四、设计报告要求

1．简述设计目的和原理。

2．按设计步骤附主要结果。

3．根据结果总结主要结论。

4.如果使用自编程序，附上源程序。

5.实验感想。

五、附录参考程序

（1）绘图程序，绘图是跟据一次实验结果绘制并非多次重复取均值后绘制

clearall

%输入信号样本个数L，FIR滤波器阶数N%

L=input（'L='）;

N=input（'N='）;

a=0.95;

%产生L个v（n）、u（n）、s（n）和x（n）,在同一座标内绘出最后100个s（n）和x（n）

w=-sqrt（3*（1-（0.95^2）））+2*sqrt（3*（1-（0.95^2）））*rand（1,L）;%w（n）的标准差

v=-sqrt（3）+2*sqrt（3）*rand（1,L）;

u=ones（1,L）;

s

（1）=1;%计算s（n）

i=2;

while（i<=L）;

s（i）=a*s（i-1）+w（i）;

i=i+1;

end;

i=1;%计算x（n）

while（i<=L）;

x（i）=s（i）+v（i）;

i=i+1;

end;

figure

（1）;%绘制s（n）（红色）,x（n）（蓝色）

k=（L-99）:

1:

L;

plot（k,s（k）,'r',k,x（k）,'k--'）;

legend（'s（n）','x（n）',0）;

title（'s（n）和x（n）'）;

xlabel（'n'）;ylabel（'输入信号'）;

%利用L个v（n）,s（n），估计Rxx和rxs,调用矩阵求逆子程序计算h1，将N个理想的h（n）和估计的h1绘与同一座标内%

fxx=zeros（1,N）;%计算fxx

for（i=1:

N）;

for（k=0:

（L-i-1））;

fxx（i）=（1/（L-i））*x（k+1）*x（k+i）+fxx（i）;

end;

end;

for（i=1:

N）;%生成Rxx矩阵

for（k=1:

N）;

Rxx（k,i）=fxx（abs（k-i）+1）;

end;

end;

rxs=zeros（1,N）;%计算rxs

for（i=1:

N）;

for（k=0:

（L-i-1））;

rxs（i）=（1/（L-i））*x（k+1）*s（k+i）+rxs（i）;

end;

end;

h1=（inv（Rxx））*rxs';

for（i=1:

N）;%绘制h（理想滤波器），h1（估计滤波器）

h（i）=0.238*（0.724）^i*u（i）;

end;

figure;

k=1:

1:

N;

plot（k,h（k）,'r',k,h1（k）,'k--'）;

title（'理想的h（n）和估算的h（n）'）;

legend（'理想的h（n）','估算的h（n）'）;

xlabel（'n'）;

%重新生成L点随机序列%

w=-sqrt（3*（1-（0.95^2）））+2*sqrt（3*（1-（0.95^2）））*rand（1,L）;

v=-sqrt（3）+2*sqrt（3）*rand（1,L）;

u=ones（1,L）;

s

（1）=1;%计算s（n）

i=2;

while（i<=L）;

s（i）=a*s（i-1）+w（i）;

i=i+1;

end;

i=1;%计算x（n）

while（i<=L）;

x（i）=s（i）+v（i）;

i=i+1;

end;

%利用公式（1－15），进行理想维纳滤波得L个SI（n）,将最后100个S（n）和SI（n）绘制于同一坐标内%

S=conv（h,x）;

SI

（1）=S

（1）;

for（i=2:

L）;

SI（i）=0.724*SI（i-1）+0.238*x（i）;

end;

figure;

k=（L-99）:

1:

L;

plot（k,s（k）,'r',k,SI（k）,'k--'）;

title（'s（n）和Si（n）'）;

legend（'s（n）','Si（n）',0）;

xlabel（'n'）;ylabel（'理想滤波输出'）;

%由h1（n）,根据公式（1-6）对x（n）进行过滤得L个SR（n）将最后100个s（n）和SR（n）绘制与同一坐标内%

SR=conv（h1,x）;

figure;

k=（L-99）:

1:

L;

plot（k,s（k）,'r',k,SR（k）,'k--'）;

title（'s（n）和Sr（n）'）;

legend（'s（n）','Sr（n）',0）;

xlabel（'n'）;ylabel（'实际滤波输出'）;

%由L个x（n）,s（n）,SI（n）,SR（n）,根据公式（1-17），统计并打印出EX2,EI2,ER2%

EX2=0;

EI2=0;

ER2=0;

for（i=1:

L）;

EX2=1/L*（（x（i）-s（i））^2）+EX2;%计算接受信号与原信号的误差

end;

for（i=1:

L）;

EI2=1/L*（（SI（i）-s（i））^2）+EI2;%计算维纳滤波器理想滤波

end;

for（i=1:

L）;%计算N点的维纳滤波器理想滤波

EI21=1/L*（（S（i）-s（i））^2）+EI2;

end;

for（i=1:

L）;%计算估计滤波器滤波

ER2=1/L*（（SR（i）-s（i））^2）+ER2;

end;

EX2

EI2

ER2

（2）检验N、L对FIR滤波器的估计精度和滤波效果的影响，为尽量保证统计意义故重复了1000次后取均值。

clearall

%输入信号样本个数L，FIR滤波器阶数N%

L=input（'L='）;

N=input（'N='）;

a=0.95;

EX=0;

EI=0;

ER=0;

for（z=1:

100）;%重复1000次

%产生L个v（n）、u（n）、s（n）和x（n）,在同一座标内绘出最后100个s（n）和x（n）%

w=-sqrt（3*（1-（0.95^2）））+2*sqrt（3*（1-（0.95^2）））*rand（1,L）;%w（n）的标准差

v=-sqrt（3）+2*sqrt（3）*rand（1,L）;

u=ones（1,L）;

s

（1）=1;%计算s（n）

i=2;

while（i<=L）;

s（i）=a*s（i-1）+w（i）;

i=i+1;

end;

i=1;%计算x（n）

while（i<=L）;

x（i）=s（i）+v（i）;

i=i+1;

end;

figure

（1）;%绘制s（n）（红色）,x（n）（蓝色）

k=（L-99）:

1:

L;

plot（k,s（k）,'r',k,x（k）,'k--'）;

legend（'s（n）','x（n）',0）;

title（'s（n）和x（n）'）;

xlabel（'n'）;ylabel（'输入信号'）;

%利用L个v（n）,s（n），估计Rxx和rxs,调用矩阵求逆子程序计算h1，将N个理想的h（n）和估计的h1绘与同一座标内%

fxx=zeros（1,N）;%计算fxx

for（i=1:

N）;

for（k=0:

（L-i-1））;

fxx（i）=（1/（L-i））*x（k+1）*x（k+i）+fxx（i）;

end;

end;

for（i=1:

N）;%生成Rxx矩阵

for（k=1:

N）;

Rxx（k,i）=fxx（abs（k-i）+1）;

end;

end;

rxs=zeros（1,N）;%计算rxs

for（i=1:

N）;

for（k=0:

（L-i-1））;

rxs（i）=（1/（L-i））*x（k+1）*s（k+i）+rxs（i）;

end;

end;

h1=（inv（Rxx））*rxs';

for（i=1:

N）;

h（i）=0.238*（0.724）^i*u（i）;

end;

figure

（2）;

k=1:

1:

N;

plot（k,h（k）,'r',k,h1（k）,'k--'）;

title（'理想的h（n）和估算的h（n）'）;

legend（'理想的h（n）','估算的h（n）'）;

xlabel（'n'）;

%重新生成L点随机序列%

w=-sqrt（3*（1-（0.95^2）））+2*sqrt（3*（1-（0.95^2）））*rand（1,L）;

v=-sqrt（3）+2*sqrt（3）*rand（1,L）;

u=ones（1,L）;

s

（1）=1;%计算s（n）

i=2;

while（i<=L）;

s（i）=a*s（i-1）+w（i）;

i=i+1;

end;

i=1;%计算x（n）

while（i<=L）;

x（i）=s（i）+v（i）;

i=i+1;

end;

%利用公式（1－15），进行理想维纳滤波得L个SI（n）,将最后100个S（n）和SI（n）绘制于同一坐标内%

S=conv（h,x）;

SI

（1）=S

（1）;

for（i=2:

L）;

SI（i）=0.724*SI（i-1）+0.238*x（i）;

end;

figure（3）;

k=（L-99）:

1:

L;

plot（k,s（k）,'r',k,SI（k）,'k--'）;

title（'s（n）和Si（n）'）;

legend（'s（n）','Si（n）',0）;

xlabel（'n'）;ylabel（'理想滤波输出'）;

%由h1（n）,根据公式（1-6）对x（n）进行过滤得L个SR（n）将最后100个s（n）和SR（n）绘制与同一坐标内%

SR=conv（h1,x）;

figure（4）;

k=（L-99）:

1:

L;

plot（k,s（k）,'r',k,SR（k）,'k--'）;

title（'s（n）和Sr（n）'）;

legend（'s（n）','Sr（n）',0）;

xlabel（'n'）;ylabel（'实际滤波输出'）;

%由L个x（n）,s（n）,SI（n）,SR（n）,根据公式（1-17），统计并打印出EX2,EI2,ER2%

EX2=0;

EI2=0;

ER2=0;

for（i=1:

L）;

EX2=1/L*（（x（i）-s（i））^2）+EX2;%计算接受信号与原信号的误差

end;

for（i=1:

L）;

EI2=1/L*（（SI（i）-s（i））^2）+EI2;%计算维纳滤波器理想滤波

end;

for（i=1:

L）;%计算估计滤波器滤波

ER2=1/L*（（SR（i）-s（i））^2）+ER2;

end;

temp1（z）=EX2;

temp2（z）=EI2;

temp3（z）=ER2;

EX=EX+EX2;

EI=EI+EI2;

ER=ER+ER2;

end;

EX=EX/1000;

EI=EI/1000;

ER=ER/1000;

EX

EI

ER