基于最小均方误差(MMSE)估计的因果维纳滤波的实现.

基于最小均方误差(MMSE)估计的因果维纳滤波的实现
一．功能简介
基于最小均方误差(MMSE)估计的因果维纳滤波的Matlab 实现，用莱文森-德宾(Levinson-Durbin)算法求解维纳-霍夫方程(Yule-wa1ker)方程，得到滤波器系数，进行维纳滤波。

二．维纳滤波简介
信号处理的实际问题，常常是要解决在噪声中提取信号的问题，因此，我们需要寻找一种所谓有最佳线性过滤特性的滤波器，这种滤波器当信号与噪声同时输入时，在输出端能将信号尽可能精确地重现出来，而噪声却受到最大抑制。

维纳(Wiener)滤波就是用来解决这样一类从噪声中提取信号问题的一种过滤(或滤波)方法。

一个线性系统，如果它的单位样本响应为h (n )，当输入一个随机信号x (n )，且
)()()(n n s n x υ+=
其中s (n )表示信号，)(n υ表示噪声，则输出y (n )为
∑-=m
m n x m h n y )()()(
我们希望x (n )通过线性系统h (n )后得到的y (n )尽量接近于s (n )，因此称y (n )为s (n )的估计值，
用)(ˆn s
表示，即 )(ˆ)(n s
n y =
维纳滤波器的输入—输出关系
如上图所示。

这个线性系统)(⋅h 称为对于s(n)的一种估计器。

如果我们以s
s ˆ与分别表示信号的真值与估计值，而用e (n )表示它们之间的误差，即
)(ˆ)()(n s
n s n e -= 显然，e (n )可能是正的，也可能是负的，并且它是一个随机变量。

因此，用它的均方值来表达误差是合理的，所谓均方误差最小即它的平方的统计平均值最小：
[][]
22)ˆ()(s
s E n e E -=最小 已知希望输出为：
1
ˆ()()()()N m y n s
n h m x n m -===-∑ 误差为：
1
ˆ()()()()()()N m e n s n s
n s n h m x n m -==-=--∑ 均方误差为：
1
2
20()(()()())N m E e n E s n h m x n m -=⎡⎤⎡⎤=--⎢⎥⎣⎦
⎣⎦
∑ 上式对() m=0,1,,N-1h m 求导得到：
1
02(()()())()0
0,1,21N opt m E s n h m x n m x n j j N -=⎡⎤---==-⎢⎥⎣⎦
∑
进一步得：
[][]
1
()()()()()0,1,1N opt m E s n x n j h m E x n m x n j j N -=-=--=-∑
从而有：
1
()()()
0,1,2,,1N xs opt xx m R j h m R j m j N -==-=-∑
于是就得到N 个线性方程：
(0)(0)(0)(1)(1)(1)(1)1(1)(0)(1)(1)(0)(1)(2)1(1)(0)(1)(1)(2)(1)(0)
xs xx xx xx xs xx xx xx xs xx xx xx j R h R h R h N R N j R h R h R h N R N j N R N h R N h R N h N R ==+++--⎧⎪==+++--⎪
⎨
⎪⎪=--=-+-++-⎩
写成矩阵形式为：
(0)
(1)(1)(0)(0)(1)(0)(2)(1)(1)(1)(2)
(0)(1)(1)xx xx xx xs xx xx xx xs xx xx xx xs R R R N R h R R R N R h R N R N R R N h N -⎡⎤⎡⎤
⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢
⎥⎢
⎥⎢⎥----⎣⎦⎣⎦⎣⎦
简化形式：xx xs R H R = 其中：H=[h(0) h(1)
h(N-1)]'是滤波器的系数
[](0),(1),(1)'xs xs xs xs R R R R N =-是互相关序列
(0)
(1)(1)(1)(0)(2)(1)(2)
(0)xx xx xx xx
xx xx xx xx xx xx R R R N R R R N R R N R N R -⎡⎤⎢⎥-⎢
⎥=⎢⎥⎢
⎥
--⎣⎦
是自相关矩阵
由上可见，设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳－霍夫（Wiener －Hopf ）方程。

另外，设计维纳滤波器要求已知信号与噪声的相关函数。

三．程序求解过程
由上述可见，本程序实现的关键是在已知输入信号的自相关函数和输入信号和理想输出信号的互相关函数的情况下，求解维纳－霍夫（Wiener －Hopf ）方程，从而得到滤波器系数，再进行维纳滤波。

求解步骤： 1. 初始化值
(0)(0)/(0)xd xx a r r = (0)(1)/(0)xd xx b r r =
2. 对于j=1,2,
,M-1，进行如下计算：
1
10()()()1(0)()()
j xd xx i j xx xx i r j r j i a i temp r r j i b i -=-=--=
--∑∑
()()2() i=0,1,
j-1a i a i temp b i =-⋅
()1a j temp =
1
01
0(1)(1)()2(0)()()
j xx xx i j xx xx i r j r i b i temp r r j i b i -=-=+-+=
--∑∑
()(1)2() i=1,
j b i b i temp b j i =--⋅-
(0)2b temp =
3．滤波器系数为：()() i=0,1,
M-1h i a i =
4．利用上面的得到的滤波器对输入信号进行维纳滤波，得到输出信号。

四．函数说明
函数使用方法：y=wienerfilter(x,Rxx,Rxd,M)
参数说明：x 是输入信号，Rxx 是输入信号的自相关向量，Rxx 是输入信号和理想信号的的互相关向量，M 是维纳滤波器的长度，输出y 是输入信号通过维纳滤波器进行维纳滤波后的输出。

具体程序见Matlab 的.m 文件。

五．程序示例
加载Matlab 中的语音数据handel ，人为地加入高斯白噪声，分别计算加入噪声后信号的自相关xx R 和加入噪声后信号和理想信号的互相关xd R ，取滤波器的长度为M=500，将以上参数代入函数中进行维纳滤波，得到输出。

程序如下：
load handel %加载语音信号 d=y; d=d*8; %增强语音信号强度 d=d';
fq=fft(d,8192); %进行傅立叶变换得到语音信号频频 subplot(3,1,1); f=Fs*(0:4095)/8192;
plot(f,abs(fq(1:4096))); %画出频谱图 title('原始语音信号的频域图形');
xlabel('频率f');
ylabel('FFT');
[m,n]=size(d);
x_noise=randn(1,n); %（0，1）分布的高斯白噪声
x=d+x_noise; %加入噪声后的语音信号
fq=fft(x,8192); %对加入噪声后的信号进行傅立叶变换，看其频谱变化
subplot(3,1,2);
plot(f,abs(fq(1:4096))); %画出加入噪声后信号的频谱图
title('加入噪声后语音信号的频域图形');
xlabel('频率f');
ylabel('FFT');
yyhxcorr=xcorr(x(1:4096)); %求取信号的信号的自相关函数
size(yyhxcorr);
A=yyhxcorr(4096:4595);
yyhdcorr=xcorr(d(1:4096),x(1:4096)); %求取信号和理想信号的互相关函数
size(yyhdcorr);
B=yyhdcorr(4096:4595);
M=500;
yyhresult=wienerfilter(x,A,B,M); %进行维纳滤波
yyhresult=yyhresult(300:8192+299);
fq=fft(yyhresult); %对维纳滤波的结果进行傅立叶变换，看其频谱变化subplot(3,1,3);
f=Fs*(0:4095)/8192;
plot(f,abs(fq(1:4096))); %画出维纳滤波后信号的频谱图
title('经过维纳滤波后语音信号的频域图形');
xlabel('频率f');
ylabel('FFT');
求出的频谱图如下所示：
由上述结果可见，经过维纳滤波后信号的噪声减弱，信噪比提高。