自适应滤波：维纳滤波器——FIR及IIR设计
作者：桂。
时间：2017-03-23 06:28:45
链接：http://www.cnblogs.com/xingshansi/p/6603263.html
【读书笔记02】
前言
仍然是西蒙.赫金的《自适应滤波器原理》第四版，距离上次看这本书已经过去半个月，要抓点紧了。本文主要包括：
1）何为维纳滤波器（Wiener Filter）;
2）Wiener滤波器的推导；
3）应用实例；
4）Wiener变体；
内容为自己的学习总结，内容多有参考他人，最后一并给出链接。
一、维纳滤波器简介
A-基本概念
对于滤波器的具体实现，都依赖两个选择：
1）Filter的impulse选择（FIR / IIR）
2）统计优化准则的选择
维纳滤波器：由数学家维纳（Rorbert Wiener）提出的一种以最小平方（统计准则）为最优准则的线性滤波器。何为线性滤波器？文中图：
B-几个问题
在具体讨论之前，先来说说自己看的时候，想到的几个问题：
问题1：维纳滤波器与自适应滤波是什么关系，与LMS呢？
个人观点：常用的hamming、taylor都是固定的滤波器，滤波器参数需要根据信号进行计算调整的，这一类滤波器都是自适应滤波器，跟自己有关嘛，自适应理所当然，维纳滤波器是利用统计参数，实际应用中可能无法得到，需要借助迭代实现，这样就成了一个新框架→自适应滤波器。参数估计呢可以利用梯度下降法迭代逼近，例如常说的LMS就是迭代逼近的一种方式，LMS本身不同于weiner-filter，但迭代结果可以作为wiener-filter的近似。
问题2：维纳滤波器是有限长（FIR-Finite Impulse Response）还是无限长滤波器（IIR-Infinite Impulse Response）？
个人观点：维纳滤波器是大的理论框架，而FIR/IIR只是实现理论的不同途径，故二者均可，下文会一一介绍。
问题3：基于维纳滤波器这个理论框架的应用有哪些？
个人观点：1）自适应是一种解决工程问题的途径，故很多自适应滤波本质都是维纳滤波（不是全部）；2）MVDR谱估计（Minimum variance distortionless response）、广义旁瓣相消GSC（Generalized Sidelobe Canceller ）等都是维纳框架下的应用。
问题4：卡尔曼滤波（Kalman）是维纳滤波的一种吗？
个人观点：应该不是，维纳滤波器是线性滤波器，而卡尔曼滤波器据说是非线性，具体区别等看到卡尔曼再回头总结。
二、维纳滤波器理论分析
A-有限长维纳滤波（FIR）
1-基本定义
下图是使用M抽头FIR滤波器的结构图：
输出为：
$\hat d\left( n \right) = \sum\limits_{k = 0}^{M - 1} {{h_k}y\left( {n - k} \right)} $
其中$h_k$为FIR滤波器系数，M为系数个数。
以下推导基于宽平稳假设，首先计算估计误差：
$e\left( n \right) = d\left( n \right) - \hat d\left( n \right) = d\left( n \right) - {{\bf{h}}^T}{\bf{y}}$
其中${{\bf{h}}^T} = \left[ {{h_0},{h_1},{h_2},...,{h_{M - 1}}} \right]$，${{\bf{y}}^T} = \left[ {y\left( n \right),y\left( {n - 1} \right),y\left( {n - 2} \right),...,y\left( {n - M + 1} \right)} \right]$是包括过去M个样本的输入向量。
Wiener Filter基于最小均方误差准则，给出均方误差定义：
其中，${\bf{r}}_{yd}^{ - 1} = E\left[ {{\bf{y}}d\left( n \right)} \right]$为输入信号与期望信号的互相关。
2-维纳滤波器求解
最小化估计误差，对其中某个抽头系数求偏导：
$\frac{{\partial J}}{{\partial {h_k}}} = 0\;\;\; \Rightarrow \;\;\; - 2E\left[ {e\left( n \right)y\left( {n - k} \right)} \right] = 0$
这就是正交定理（Orthogonality Principle）：估计误差$e(n)$需要正交与输入信号$y(n)$。这也容易理解，在输入信号$y(n)$张成的子空间中，更加高维的信息无法被表达，故成了误差。如$（x_1,y_1,z_1）$用$x 、 y$两个单位向量表达，最小均方误差时$z_1$就成了估计误差。
不失一般性，将偏导写成向量/矩阵形式：
$\frac{{\partial J}}{{\partial {\bf{h}}}} = 0\;\;\; \Rightarrow \;\;\; - 2{\bf{r}}_{yd}^{ - 1} + 2{{\bf{h}}^T}{{\bf{R}}_{yy}} = 0$
上式${{\bf{h}}^T}{{\bf{R}}_{yy}} = {\bf{r}}_{yd}^{ - 1}$就是Wiener Hopf方程。
得到Wiener Filter最优解${{{\bf{h}}_{opt}}}$:
${{\bf{h}}_{opt}}{\rm{ = }}{\bf{R}}_{_{yy}}^{{\rm{ - 1}}}{\bf{r}}_{yd}^{ - 1}$
上式就是Wiener-Hopf的解，也就是对应Wiener Filter的解。求解需要矩阵求逆，又相关矩阵为对称且为Toeplitz形式，故可借助数学手段对R高效求逆——Levinson-Durbin算法。
因为在时域分析，此时FIR的解也叫时域维纳滤波器。
B-无限长维纳滤波（IIR）
1-基本定义
滤波器为无限长时，$\hat d\left( n \right) = \sum\limits_{k = 0}^{M - 1} {{h_k}y\left( {n - k} \right)}$改写为：
$\hat d\left( n \right) = \sum\limits_{k = - \infty }^\infty {{h_k}y\left( {n - k} \right)} = h\left( n \right) * y\left( n \right)$
$*$表示卷积。易证：时域有限长对应频域无限长，时域无限长对应频域有限长，因此对于IIR情形，更希望在频域进行分析。
2-维纳滤波器求解
计算均方误差：
其中${E\left( {{\omega _k}} \right)}$为$e\left( n \right)$的频域变换，${P_{yd}}\left( {{\omega _k}} \right) = E\left[ {Y\left( {{\omega _k}} \right){D^*}\left( {{\omega _k}} \right)} \right]$, ${P_{yy}}\left( {{\omega _k}} \right) = E\left[ {{{\left| {Y\left( {{\omega _k}} \right)} \right|}^2}} \right]$。
针对J求解复导数：
$\frac{{\partial J}}{{\partial H\left( {{\omega _k}} \right)}} = 0\;\;\; \Rightarrow \;\;\;{H^*}\left( {{\omega _k}} \right){P_{yy}}\left( {{\omega _k}} \right) - {P_{yd}}\left( {{\omega _k}} \right) = 0$
得到频域维纳滤波最优解：
${H_{opt}}\left( {{\omega _k}} \right) = \frac{{{P_{dy}}\left( {{\omega _k}} \right)}}{{{P_{yy}}\left( {{\omega _k}} \right)}}$
因为在频域分析，此时IIR的解也叫频域维纳滤波器。
三、应用实例
已知：
含有噪声的正弦波：$y(n) = s(n) + w(n) = \sin (2\pi fn + \theta ) + w(n)$.
其中$f = 0.2$为归一化频率[-1/2, 1/2]，$\theta$为正弦波相位，服从[0,2$\pi$]的均匀分布，$w(n)$为具有零均值和方差$\sigma^2 = 2$的高斯白噪声。
求：
时域以及频域维纳滤波器。假设滤波器为时域滤波器时$M=2$.
A-对于时域维纳滤波器：
首先求解相关矩阵：
$x(n)$为广义平稳随机过程，可以计算其自相关函数：
${r_{xx}}\left( m \right) = \cos (2\pi fn)$
利用上面推导的Wiener-Hopf最优解公式：
当然也可以求解准则函数$J$,求极值点：
B-对于频域维纳滤波器：
白噪声与信号不相关，直接调用上文推导的公式：
$H\left( {{\omega _k}} \right) = \frac{{{P_{xy}}\left( {{\omega _k}} \right)}}{{{P_{yy}}\left( {{\omega _k}} \right)}} = \frac{{{P_{xx}}\left( {{\omega _k}} \right)}}{{{P_{xx}}\left( {{\omega _k}} \right) + {P_{nn}}\left( {{\omega _k}} \right)}} = \frac{{{P_{xx}}\left( {{\omega _k}} \right)}}{{{P_{xx}}\left( {{\omega _k}} \right) + 2}}$
时域相关函数与频域功率谱互为傅里叶变换，而时域相关函数在求时域维纳滤波时已写出，得出功率谱公式：
当$L -> \infty $，将$P_{xx}$代入上式，即使频域维纳滤波器的解。
C-闲话时域、频域维纳滤波器
时域求解维纳滤波器，但时域对应的都有频域的变换结果，如$L = 2$的频域解就是时域维纳解的傅里叶变换。画出$L$不同取值对应的频域幅度响应：
可以观察到：$L$趋近∞时，频域响应接近$\sigma(.)$冲激响应，这与理论：相关函数为时域余弦对应频域冲激响应是一致的。
这个例子只是理想情况，理一理求解的思路，事实上认为信号的自相关为已知，这是不符合实际的。实际应用中如何近似求解呢？给出一个简单例子。
D-基于频域维纳滤波的语音增强
还是利用上面的模型:
$y(n) = x(n) + w(n)$
这里$y(n)$是麦克风接收的带噪语音，$x(n)$是干净语音信号，$w(n)$为白噪声。显然相关函数我们无法得知。
利用一种近似的处理思路：利用前面几个分帧不带语音，估计噪声，从而得到噪声的功率谱近似，利用带噪语音功率减去噪声功率，得到
$H\left( {{\omega _k}} \right) = \frac{{{P_{xy}}\left( {{\omega _k}} \right)}}{{{P_{yy}}\left( {{\omega _k}} \right)}} = \frac{{{P_{xx}}\left( {{\omega _k}} \right)}}{{{P_{xx}}\left( {{\omega _k}} \right) + {P_{nn}}\left( {{\omega _k}} \right)}} = \frac{{{P_{yy}}\left( {{\omega _k}} \right) - {P_{nn}}\left( {{\omega _k}} \right)}}{{{P_{yy}}\left( {{\omega _k}} \right)}}$
利用估计出的维纳滤波器，即可实现信号的频域滤波。这里只是想到的一个实际例子，至于参数估计、迭代方式则是百花齐放了。
附上主要代码：
nw = 512; ni = 64;
NIS = 100;
Y = fft(enframe(y,hamming(nw),ni)');
Yab = abs(Y);
Nest = mean(Yab(:,1:NIS),2);
Yest = zeros(size(Y));
for i = 1:size(Y,2)
Yest(:,i) = Yab(:,i).*((Yab(:,i)-Nest)./Yab(:,i));
end
Ye = Yest.*exp(1j*angle(Y));
result = zeros(1,length(y));%estimation
for i =1:size(Y,2);
pos = ((i-1)*ni+1):((i-1)*ni+nw);
result(pos) = result(pos)+real(ifft(Ye(:,i))).';
end
result = result/max(abs(result));
上面思路处理结果：
可以看出维纳降噪多少还是有些效果的，$H\left( {{\omega _k}} \right) = \frac{{{P_{xy}}\left( {{\omega _k}} \right)}}{{{P_{yy}}\left( {{\omega _k}} \right)}} = \frac{{{P_{xx}}\left( {{\omega _k}} \right)}}{{{P_{xx}}\left( {{\omega _k}} \right) + {P_{nn}}\left( {{\omega _k}} \right)}} = \frac{1}{{1 + 1/SNR}}$可以看出SNR越小，维纳滤波器衰减越大。
四、Wiener Filter变体
A-平方根维纳滤波
使用维纳滤波器的平方根，则滤波器在频域的滤波结果为：
$\hat X\left( {{\omega _k}} \right) = \sqrt {H\left( {{\omega _k}} \right)} Y\left( {{\omega _k}} \right)$
仍然基于噪声与信号不相关的假设，分析滤波后信号的功率谱：
$E{\left| {\hat X\left( {{\omega _k}} \right)} \right|^2} = {P_{\hat x\hat x}} = H\left( {{\omega _k}} \right)E{\left| {Y\left( {{\omega _k}} \right)} \right|^2} = H\left( {{\omega _k}} \right){P_{yy}}\left( {{\omega _k}} \right) = {P_{xx}}\left( {{\omega _k}} \right)$
可见采用平方根维纳滤波，滤波器输出信号的功率谱与纯净信号的功率谱相等。
B-参变维纳滤波器
以频域Wiener Filter为例：
$H\left( {{\omega _k}} \right) = \frac{{{P_{xy}}\left( {{\omega _k}} \right)}}{{{P_{yy}}\left( {{\omega _k}} \right)}} = \frac{{{P_{xx}}\left( {{\omega _k}} \right)}}{{{P_{xx}}\left( {{\omega _k}} \right) + {P_{nn}}\left( {{\omega _k}} \right)}}$
很自然地，可以将其扩展为广义Wiener Filter形式：
$H\left( {{\omega _k}} \right) = {\left( {\frac{{{P_{xx}}\left( {{\omega _k}} \right)}}{{{P_{xx}}\left( {{\omega _k}} \right) + \alpha {P_{nn}}\left( {{\omega _k}} \right)}}} \right)^\beta }$
这样通过调节（$\alpha$,$\beta $ ）即可对滤波器进行微调。
参考：
Philipos C.Loizou《speech enhancement theory and practice》.
Simon Haykin 《Adaptive Filter Theory Fourth Edition》.
