作者：桂。

時間：2017-03-23  06:28:45

 【讀書筆記02】

前言

仍然是西蒙.赫金的《自適應濾波器原理》第四版，距離上次看這本書已經過去半個月，要抓點緊了。本文主要包括：

　　1）何為維納濾波器（Wiener Filter）;

　　2）Wiener濾波器的推導；

　　3）應用例項；

　　4）Wiener變體；

內容為自己的學習總結，內容多有參考他人，最後一併給出連結。

一、維納濾波器簡介

　　A-基本概念

對於濾波器的具體實現，都依賴兩個選擇：

1）Filter的impulse選擇（FIR / IIR）

2）統計優化準則的選擇

維納濾波器

　　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》.