最小均方(LMS)自适应滤波算法

1 LMS函数实现

% ----------------------------LMS(Least Mean Squre)算法------------------------------ %
% parm:
%     xn   输入的信号序列      (列向量)
%     dn   所期望的响应序列    (列向量)
%     M    滤波器的阶数        (标量)
%     mu   收敛因子(步长)      (标量)     要求大于0,小于xn的相关矩阵最大特征值的倒数    
%     itr  迭代次数            (标量)     默认为xn的长度,M<itr<length(xn)
% ---------------------------------------------------------------------------------- %
% retval:
%     W    滤波器的权值矩阵     (矩阵)
%          大小为M x itr,
%     en   误差序列(itr x 1)    (列向量)  
%     yn   实际输出序列         (列向量)
% 参数个数必须为4个或5个
% ---------------------------------------------------------------------------------- %
function [yn,W,en]=LMS(xn,dn,M,mu,itr)
%% 检查输入参数个数 
if nargin == 4                 % itr参数缺省时，将itr设为输入信号的长度
    itr = length(xn);
elseif nargin == 5             % 无参数缺省时
	if itr > length(xn) | itr < M  % 错误判定
        error('迭代次数过大或过小!');
    end
else                           % 错误判定
    error('请检查输入参数的个数!');
end

%% 创建，初始化参数
en = zeros(itr,1);             % 误差序列,en(k)表示第k次迭代时预期输出与实际输入的误差
W  = zeros(M,itr);             % 每一行代表一个加权参量,每一列代表每次迭代,初始为0

%% 迭代计算
for k = M:itr                  % 第k次迭代(从M次迭代到itr次)
    x = xn(k:-1:k-M+1);        % 提取长度为M的输入信号向量（注意MATLAB中索引是从1开始的，所以-1表示倒序），k:-1:k-M+1 表示从当前位置 k 开始向前取M个值
    y = W(:,k-1).' * x;        % 计算出当前迭代滤波器的输出信号，W(:,k-1)表示取前一次迭代的滤波器权值向量，.'表示转置，* 表示矩阵乘法
    en(k) = dn(k) - y ;        % 第k次迭代的误差（期望信号与实际输出信号间的误差）
    % 滤波器权值计算的迭代式
    W(:,k) = W(:,k-1) + 2*mu*en(k)*x; % 使用LMS算法更新滤波器权值
end

%% 求最优时滤波器的输出序列
yn = inf * ones(size(xn));
for k = M:length(xn)
    x = xn(k:-1:k-M+1);
    yn(k) = W(:,end).' * x;
end

整个函数实现主要围绕的就是两个参数，一个就是误差序列 $e n$ ，另一个就是滤波器的权值矩阵 $M$ ，然后做循环获得最好的滤波器权值矩阵，然后就用这个最优解求出最优的输出信号

核心算法就是如何更新滤波器权值

en(k) = dn(k) - y;
W(:,k) = W(:,k-1) + 2*mu*en(k)*x; % 使用LMS算法更新滤波器权值

这个算法感觉很像 $P I D 控制算法$ 和 $最短路径的搜索算法$ （不停做循环操作，计算 $error$ 参数，将 $error$ 代入算法中进行修正，一步步减小误差，进而不断逼近期望信号 $d n$ ）

2 用LMS自适应算法对语语音降噪测试

close all;
clear all; 
clc; 
 
[s, fs] = audioread('C5_1_y.wav');          % 读入数据文件

%% 对上面获取的语言信号，其实就是无噪信号进行预处理
s=s-mean(s);                                % 消除直流分量（原信号-原信号经过均值后的信号）
s=s/max(abs(s));                            % 幅值归一化处理

N=length(s);                                % 语音长度
time=(0:N-1)/fs;                            % 设置时间刻度


%% 对纯净信号添加噪声是之无噪信号变成含噪信号
SNR=5;                                      % 设置信噪比
r1   = awgn(s,SNR,'measured','db');		    % 添加加性白高斯噪声（Additive White Gaussian Noise，AWGN）后得到含噪信号r1
itr  = length(r1);							% 有没有这句不重要，LMS函数实现对缺省情况拟以代码处理这类情况  
snr1 = SNR_Calc(s,r1);                      % 计算初始信噪比

%% 设置降噪过程需要的固定参数
M=64;                                       % 设置Ｍ和mu
mu=0.001;

[y,W,e] = LMS(r1,s,M,mu,itr);

output = e/max(abs(e));                 	% LMS滤波输出
snr2=SNR_Calc(s,output);                    % 计算滤波后的信噪比

snr=snr2-snr1;								% 得到滤波前后的信噪比差值

SN1=snr1; SN2=snr2; SN3=snr;

fprintf('snr1=%5.4f   snr2=%5.4f    snr=%5.4f\n',snr1,snr2,snr);

% 作图
subplot 311; plot(time,s,'k'); ylabel('幅值') 
ylim([-1 1]); title('原始语音信号');
subplot 312; plot(time,r1,'k'); ylabel('幅值') 
ylim([-1 1]); title('带噪语音信号');
subplot 313; plot(time,output,'k'); 
ylim([-1 1]); title('LMS滤波输出语音信号');