作者：紫曜花
時間：2018-11-30
前言
姚天任、孫洪的《現代數字信號處理》第三章自適應濾波中關於LMS算法的學習，全文包括：
1. 自適應濾波器簡介
2. 自適應干擾抵消原理
3. 自適應濾波原理
4. 最小均方（LMS）算法
5. Matlab實現
內容為自己讀書記錄，本人知識有限，若有錯誤之處，還請各位指出！
一、自適應濾波器簡介
自適應濾波器由參數可調的數字濾波器和自適應算法兩部分組成。如圖所示。
輸入信號x(n) 通過參數可調數字濾波器后產生輸出信號 y(n)，將其與期望信號d(n)進行比較，形成誤差信號e(n)， 通過自適應算法對濾波器參數進行調整，最終使 e(n)的均方值最小。自適應濾波可以利用前一時刻已得的濾波器參數的結果，自動調節當前時刻的濾波器參數，以適應信號和噪聲未知的或隨時間變化的統計特性，從而實現最優濾波。自適應濾波器實質上就是一種能調節自身傳輸特性以達到最優的維納濾波器。自適應濾波器不需要關於輸入信號的先驗知識，計算量小，特別適用於實時處理。維納濾波器參數是固定的，適合於平穩隨機信號。卡爾曼濾波器參數是時變的，適合於非平穩隨機信號。然而，只有在信號和噪聲的統計特性先驗已知的情況下，這兩種濾波技術才能獲得最優濾波。在實際應用中，常常無法得到信號和噪聲統計特性的先驗知識。在這種情況下，自適應濾波技術能夠獲得極佳的濾波性能，因而具有很好的應用價值。
常用的自適應濾波技術有：最小均方（LMS）自適應濾波器、遞推最小二乘（RLS）濾波器格型濾波器和無限沖激響應（IIR）濾波器等。這些自適應濾波技術的應用又包括：自適應噪聲抵消、自適應譜線增強和陷波等。
二、 自適應干擾抵消原理
如圖是自適應干擾抵消器的基本結構。其中期望響應d(n) 是信號與噪聲之和，即
，自適應處理器的輸入是與 相關的另一個噪聲 。當 
與 
不相關時，利用噪聲的相關性,自適應濾波器將調整自己的參數，以力圖使
成為 
的最佳估計 。這樣， 
將逼近信號
。噪聲 
就得到一定程度的抵消。但若參考通道除檢測到噪聲 
外，還收到信號分量，則自適應濾波器的輸出中將包含信號分量，從而使噪聲對消效果變壞。因此，為獲得良好的噪聲對消性能，應使參考通道檢測到的信號盡可能小，在信號不可測的噪聲環境拾取參考輸入信號。
三、自適應濾波原理
單輸入情況：輸入信號矢量 
輸出 
（3.1）
自適應線性組合器的L+1個權系數構成一個權系數矢量，稱為權矢量，用
表示，
即
則（3.1）可表示為
（3.2）
參考響應與輸出響應之差稱為誤差信號，用
表示，即
自適應線性組合器按照誤差信號均方值（或平均功率）最小的准則，即
由式（3.3）和（3.4）可得均方誤差表示式為
在
和
都是平穩隨機信號的情況下，有如下定義：
輸入信號
的自相關矩陣
期望信號與輸入信號的互相關矩陣
則均方誤差的簡單表示形式
從該式可看出，在輸入信號和參考響應都是平穩隨機信號的前提下，均方誤差是權矢量的各分量的二次函數。的函數圖形是L+2維空間中一個中間下凹的超拋物面，有唯一的最低點,該曲面稱為均方誤差性能曲面，簡稱性能曲面。
均方誤差性能曲面的梯度
令梯度
等於零，可求得最小均方誤差對應的最佳權矢量或維納解
。
解得
四、最小均方(LMS)算法
最陡下降法是沿性能曲面最陡方向向下搜索曲面的最低點。曲面的最陡下降方向是曲面的負梯度方向。這是一個迭代搜索過程。最陡下降法迭代計算權矢量的公式為
最陡下降法每次迭代都需要知道性能曲面上某點的梯度值，而實際上梯度值只能根據觀測數據進行估計。LMS算法是一種很有用且很簡單的估計梯度的方法。只需要每次迭代運算時知道輸入響應和參考信號。
LMS算法最核心的思想是用平方誤差代替均方誤差。
實際上，
只是單個平方誤差序列的梯度，而
則是多個平方誤差序列統計平均的梯度，所以LMS算法就是用前者作為后者的近視。由(4.1)可得
該式是LMS算法的基本關系式。該式說明，LMS算法實際上是在每次迭代中使用粗略的梯度估計值
來代替精確值
。不難預計，權系數的調整路徑不可能准確地沿着理想的最陡下降的路徑，因而權系數的調整過程是“有噪”的，或者說
不再是確定性函數而變成了隨機變量。
下一刻權矢量
等於當前權矢量
加上一個修正量，該修正量等於誤差信號
的加權值，加權系數為
，它正比於當前的輸入信號。
五、Matlab實現
程序借用自 LMS算法自適應濾波器
1.LMS算法實現
% 迭代計算
for k = M:itr % 第k次迭代
x = xn(k:-1:k-M+1); % 濾波器M個抽頭的輸入
y = W(:,k-1).' * x; % 濾波器的輸出
en(k) = dn(k) - y ; % 第k次迭代的誤差
% 濾波器權值計算的迭代式
W(:,k) = W(:,k-1) + 2*mu*en(k)*x;
end
LMSfilter.m
% 輸入參數:
% xn 輸入的信號序列 (列向量)
% dn 所期望的響應序列 (列向量)
% M 濾波器的階數 (標量)
% mu 收斂因子(步長) (標量) 要求大於0,小於xn的相關矩陣最大特征值的倒數
% 輸出參數:
% W 濾波器的權值矩陣 (矩陣)
% 大小為M x itr,
% en 誤差序列(itr x 1) (列向量)
% yn 實際輸出序列 (列向量)
function [yn,W,en]=LMSfilter(xn,dn,M,mu)
itr = length(xn);
en = zeros(itr,1); % 誤差序列,en(k)表示第k次迭代時預期輸出與實際輸入的誤差
W = zeros(M,itr); % 每一行代表一個加權參量,每一列代表-次迭代,初始為0
% 迭代計算
for k = M:itr % 第k次迭代
x = xn(k:-1:k-M+1); % 濾波器M個抽頭的輸入
y = W(:,k-1).' * x; % 濾波器的輸出
en(k) = dn(k) - y ; % 第k次迭代的誤差
% 濾波器權值計算的迭代式
W(:,k) = W(:,k-1) + 2*mu*en(k)*x;
end
% 求最優時濾波器的輸出序列 r如果沒有yn返回參數可以不要下面的
yn = inf * ones(size(xn)); % inf 是無窮大的意思
for k = M:length(xn)
x = xn(k:-1:k-M+1);
yn(k) = W(:,end).'* x;%用最后得到的最佳估計得到輸出
end
2. 主函數main,m設計
產生信號源，可以直接讀取音頻文件。
產生相關噪聲，一個作為濾波輸入，一個產生期望信號。
產生期望信號，信號源於噪聲疊加產生含噪信號作為期望信號。
LMS濾波算法，調用LMSfilter.m函數。
繪制圖形進行觀察對比。
clc;
clear all;
close all;
%% 產生信號源
[X,Fs] = audioread('voicefile.wav');
s = X(:,1); %取出雙通道中其中一個通道作為信號源s
audiowrite('原始音頻.wav',s,Fs); %創建原始音頻.wav
n = length(s);
t=(0:n-1);
figure(1);
subplot(4,1,1);
plot(t,s);grid;ylim([-2 2]);
ylabel('幅度');
xlabel('時間');
title('原始音頻信號');
%% 產生均值為0方差為0.1的噪聲信號
v = sqrt(0.1)*randn(n,1);
%% 產生AR模型的噪聲
ar=[1,1/2]; %AR模型
v_ar=filter(1,ar,v);
% subplot(4,1,2);
% plot(t,v_ar);grid;
% ylabel('幅度');
% xlabel('時間');
% title('AR模型噪聲信號');
%% 產生MA模型的噪聲 是AR模型的相關噪聲
ma=[1,-0.8,0.4,-0.2]; %MA模型
v_ma=filter(ma,1,v);
subplot(4,1,2);
plot(t,v_ma);grid;ylim([-2 2]);
ylabel('幅度');
xlabel('時間');
title('相關噪聲信號');
%% 產生期望信號
dn = s + v_ar;
audiowrite('含噪音頻.wav',dn,Fs); %創建含噪音頻
subplot(4,1,3);
plot(t,dn);grid;ylim([-2 2]);
ylabel('幅度');
xlabel('時間');
title('含噪音頻信號');
%% LMS濾波算法
M = 50; %濾波器階數M
mu = 0.0008; %濾波器的步長
[ylms,W,elms] =LMSfilter(v_ma,dn,M,mu);
%% 繪制去噪后的語音信號
subplot(4,1,4);
plot(t,elms);grid;ylim([-2 2]);
ylabel('幅度');
xlabel('時間');
title('去噪后的音頻信號');
audiowrite('去噪音頻.wav',elms,Fs);%保存去除噪聲的音頻
%%
e = s-elms;%剩余噪聲
figure(2);
subplot(2,1,1);
plot(t,e);grid;
ylabel('幅度');
xlabel('時間');
title('剩余噪聲');
%% 一小段三個信號比較
subplot(2,1,2);
t=(100000:100500);
plot(t,elms(100000:100500,1 ),'r',t,e(100000:100500,1),'g',t,s(100000:100500,1),'b');
axis([100000,100500,-1,1]);
ylabel('幅度');
xlabel('時間');
legend('去噪后的語音信號','剩余噪聲','原始音頻');
title('一小段三個信號比較');
3.結果分析
步長 $\mu =0.0001$時
步長$\mu =0.0004$時
步長$\mu =0.0012$時
從圖中可以看出雖然真實的音頻信號中混雜了很強的噪聲，甚至噪聲淹沒了真實的信號，但是通過我們的LMS自適應濾波器后，可以很好的恢復出真實信號。所以自適應濾波器具有良好的去噪性能。
通過改變步長 $\mu$，可以發現隨着步長增加，收斂時間t逐漸減少，音頻信號在$ \mu$ 很小的時候，因為收斂慢，混有噪聲；但是當$\mu$ 增加是，音頻信號越來越清晰。但是當步長超過某一值時，音頻中又開始很有噪聲產生，再接着增加步長$\mu$，甚至會把音頻信號也濾去掉。所以若 $\mu$ 值取的過小，收斂速度就會過於緩慢，當 $\mu$ 取的過大時，又會造成系統收斂的不穩定，導致發散。
參考：
姚天任、孫洪《現代數字信號處理》