作者：桂。
時間：2017-04-02 08:08:31
鏈接：http://www.cnblogs.com/xingshansi/p/6658203.html
聲明：歡迎被轉載，不過記得注明出處哦~
【讀書筆記08】
前言
西蒙.赫金的《自適應濾波器原理》第四版第五、六章：最小均方自適應濾波器（LMS，Least Mean Square）以及歸一化最小均方自適應濾波器（NLMS，Normalized Least Mean Square）。全文包括：
1）LMS與維納濾波器（Wiener Filter）的區別；
2）LMS原理及推導；
3）NLMS推導；
4）應用實例；
內容為自己的讀書記錄，其中錯誤之處，還請各位幫忙指出！
一、LMS與維納濾波器（Wiener Filter）的區別
- 這里介紹的LMS/NLMS，通常逐點處理，對應思路是：隨機梯度下降；
- 對於Wiener Filter，給定准則函數J，隨機/批量梯度都可以得出最優解；
- LMS雖然基於梯度下降，但准則僅僅是統計意義且通常引入誤差，可以定義為$J_0$，簡而言之$J$通常不等於$J_0$，得出的最優解$w_o$自然也通常不等於維納最優解；
- 分析LMS通常會分析穩定性，穩定性是基於Wiener解，之前已給出分析。但LMS是Wiener解的近似，所以：迭代步長的穩定性，嚴格適用於Wiener解，對於LMS只是一種近似參考，並沒有充分的理論依據。
下文的分析仍然隨機梯度下降的思路進行。
二、LMS原理及推導
LMS是時間換空間的應用，如果迭代步長過大，仍然有不收斂的問題；如果迭代步長過小，對於不平穩信號，還沒有實現尋優就又引入了新的誤差，屋漏偏逢連夜雨！所以LMS系統是脆弱的，信號盡量平穩、哪怕短時平穩也湊合呢。
給出框圖：
關於隨機梯度下降，可以參考之前的文章。這里直接給出定義式：
利用梯度下降：
$- \nabla J = {\bf{x}}{\left( {{{\bf{w}}^T}{\bf{x}} - {d}} \right)^T}$
給出LMS算法步驟：
1）給定$\bf{w}(0)$，且$1<\mu<1/\lambda_{max}$；
2）計算輸出值：$y\left( k \right) = {\bf{w}}{\left( k \right)^T}{\bf{x}}\left( k \right)$;
3）計算估計誤差：$e\left( k \right) = d\left( k \right) - y\left( k \right)$;
4）權重更新：${\bf{w}}\left( {k + 1} \right) = {\bf{w}}\left( k \right) + \mu e\left( k \right){\bf{x}}\left( k \right)$
三、NLMS推導
看到Normalized，與之聯系的通常是約束條件，看到約束不免想起拉格朗日乘子。思路有了，現在開始分析：
假設${\bf{w}}\left( k \right) \Rightarrow {\bf{w}}\left( {k + 1} \right)$得到最優權重，即：
$d\left( k \right) = {\bf{w}}\left( {k + 1} \right){\bf{x}}\left( k \right)$
我們希望在得到期望權重的附近，迭代不要過大以免錯過最優值：
寫出准則函數：
利用之前文章提到的拉格朗日乘子法：
這里僅僅分析基於歐式距離$p = 2$的情形，其它范數類似。求解得出:
通常為了防止分母為零迭代方程需要修正，而修正后步長存在偏差，故添加調節因子$\mu$：
給出NLMS算法步驟：
1）給定$\bf{w}(0)$；
2）計算輸出值：$y\left( k \right) = {\bf{w}}{\left( k \right)^T}{\bf{x}}\left( k \right)$;
3）計算估計誤差：$e\left( k \right) = d\left( k \right) - y\left( k \right)$;
4）權重更新：${\bf{w}}\left( {k + 1} \right) = {\bf{w}}\left( k \right) + \frac{\mu }{{\alpha + {{\left| {{\bf{x}}\left( k \right)} \right|}^2}}}{\bf{x}}\left( k \right){e^*}\left( k \right)$
四、應用實例
A-自適應噪聲濾波
這個場景可以簡化為：一個房間兩個麥克風，一個放在遠處采集房間噪聲，一個放在說話人附近采集帶噪語音，認為兩個音頻文件的噪聲相似。
這里噪聲直接用白噪聲，對應實際場景可以認為是采集的噪聲數據，給出主要代碼：
[s, fs, bits] = wavread(filename); s=s-mean(s); s=s/max(abs(s)); N=length(s); time=(0:N-1)/fs; %%生成帶噪信號 clean=s'; ref_noise=0.1*randn(1,length(s)); mixed = clean+ref_noise %NLMS mu=0.05;M=32;espon=1e-4; % [en,wn,yn]=lmsFunc(mu,M,ref_noise,mixed);% [en,wn,yn]=nlmsFunc(mu,M,ref_noise,mixed,espon);
LMS代碼：
function [e,w,ee]=lmsFunc(mu,M,u,d)
% Normalized LMS
% Call:
% [e,w]=nlms(mu,M,u,d,a);
%
% Input arguments:
% mu = step size, dim 1x1
% M = filter length, dim 1x1
% u = input signal, dim Nx1
% d = desired signal, dim Nx1
% a = constant, dim 1x1
%
% Output arguments:
% e = estimation error, dim Nx1
% w = final filter coefficients, dim Mx1
%intial value 0
w=zeros(M,1); %This is a vertical column
%input signal length
N=length(u);
%make sure that u and d are colon vectors
u=u(:);
d=d(:);
%NLMS
ee=zeros(1,N);
for n=M:N %Start at M (Filter Length) and Loop to N (Length of Sample)
uvec=u(n:-1:n-M+1); %Array, start at n, decrement to n-m+1
e(n)=d(n)-w'*uvec;
w=w+2*mu*uvec*e(n);
% y(n) = w'*uvec; %In ALE, this will be the narrowband noise.
end
NLMS代碼：
function [e,w,ee]=nlmsFunc(mu,M,u,d,a)
% Normalized LMS
% Call:
% [e,w]=nlms(mu,M,u,d,a);
%
% Input arguments:
% mu = step size, dim 1x1
% M = filter length, dim 1x1
% u = input signal, dim Nx1
% d = desired signal, dim Nx1
% a = constant, dim 1x1
%
% Output arguments:
% e = estimation error, dim Nx1
% w = final filter coefficients, dim Mx1
%intial value 0
w=zeros(M,1); %This is a vertical column
%input signal length
N=length(u);
%make sure that u and d are colon vectors
u=u(:);
d=d(:);
%NLMS
ee=zeros(1,N);
for n=M:N %Start at M (Filter Length) and Loop to N (Length of Sample)
uvec=u(n:-1:n-M+1); %Array, start at n, decrement to n-m+1
e(n)=d(n)-w'*uvec;
w=w+mu/(a+uvec'*uvec)*uvec*e(n);
% y(n) = w'*uvec; %In ALE, this will be the narrowband noise.
end
對應結果圖：
可以看出LMS/NLMS在最開始都有一個自適應的過程。
NLMS基於信號$x$的能量實現變步長，信號大步長小，信號小則步長大：目標信號明顯，則迭代細致，不明顯，則一帶而過，呵呵，跟平時看書還挺像，聰明的孩子。
再來看一組信號：
這里在中間令噪聲突變，可以看到濾波器又需要重新自適應，因此對於短時平穩LMS勉強使用，如果不斷變呢？非平穩LMS自然無效了，這個時候就需要Kalman Filter來搭把手。
B-工頻噪聲濾波
現在有一個音頻信號，分析頻譜：
可以看到信號帶有明顯的$50Hz$噪聲，我們知道$50Hz$的正弦與余弦可以組合成任意相位的$50Hz$頻率信號，基於這個思路，進行自適應濾波：
給出主要的代碼：
x1=cos(2*pi*50*time); x2=sin(2*pi*50*time); w1=0.1; w2=0.1; e=zeros(1, N); y=zeros(1, N); mu=0.05; for i=1: N y(i)=w1 * x1(i)+ w2 * x2(i); e(i) =x(i)-y(i); w1=w1+mu * e(i) * x1(i); w2=w2+mu * e(i) * x2(i); end
結果圖可以看出，工頻50Hz濾除：
基於LMS的應用還有很多，不一一說啦。
參考：
- Simon Haykin 《Adaptive Filter Theory Fourth Edition》.
- 宋知用：《MATLAB在語音信號分析和合成中的應用》.