作者：桂。
時間：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》.