自适应滤波器减噪

0. 前言知识

什么是自适应滤波？

自适应滤波：根据前一时刻的序列，自动的调节这一时刻的滤波器参数。

实质：一种能调节自身传输特性以达到最优的维纳滤波器。

几种自适应滤波器的比较

自适应：无需信号的先验知识，适用于实时处理。
维纳滤波：参数固定，适用于平稳随机信号。
卡尔曼滤波器：参数时变，适合于非平稳随机信号。

常用滤波器算法

LMS(最小均方)自适应滤波器
RLS(递推最小二乘)滤波器
IIR(无限冲激相应)滤波器

应用分类：

自适应噪声抵消
自适应谱线增强
陷波

1. LMS滤波

LMS原理及推导

$X(n) = [x(n),x(n-1),\dots,x(n-L)]$ 是一截输入信号序列， $d(n)$ 为期望输出信号，两者之差定义为误差序列：

$e(n) = d(n) - \sum_{i=1}^L w_ix(n-i)$

LMS中滤波器的权系数则为： $w_i$

注：有的博客或者书中权重系数都是 $w_i(n)$ 代表权重是动态变化的。

用矩阵表达：

$e(n) = d(n) - W^TX(n)$

误差二次方为：

$e^2(n) = d^2(n) - 2d(n)X^T(n)W + W^TX(n)X^T(n)W$

自适应线性组合器按照误差信号均方值最小的准则，即

$E[e^2(n)] = E[d^2(n)] -2E[d(n)X^T(n)]W+W^TE[X(n)X^T(n)]W$

定义自相关矩阵R和互相关矩阵P

$\begin{aligned} R = & E[X(n)X^T(n)] \\ P = & E[d(n)X^T(n)] \end{aligned}$

则有：

$E[e^2(n)] = E[d^2(n)] + W^TRW - 2P^TW$

在 $d(n)$ 和 $X(n)$ 都是平稳随机信号的情况下，均方误差是权矢量( $W$ ,由 $L+1$ 个权系数 $w_i$ 构成)的各分量的二次函数。该曲面是一个 $L+2$ 维空间的上凸下凹的超抛物面。（这还需待补充，画张图）

这里还需要补充随机信号的时间平均与集平均的概念

集平均就是指多次同一个随机信号多次测试取平均。根据滤波器的宽度，这里的期望就是多次平移信号取均值。

求一阶导数：

$\nabla(n) =\frac{\partial{E[e^2(n)]}}{\partial{W}} = 2RW - 2P$

令梯度等于0，可以求得最小均方误差对应的最佳权矢量 $W$ 或维纳解 $W^*$ 。

$W^* = R^{-1}P$

最小均方(LMS)算法

若直接按照上面的公式去计算，需要精确地知道输入信号和期望信号的先验统计知识 $R$ ，这就无法做到实时的滤波效果，为了避免求逆操作，这里采用迭代的方式搜索曲面的最低点。曲面最陡下降方向是负梯度方向。最陡下降法迭代计算权矢量公式为：

$W(n+1) = W(n) - \mu\nabla(n)$

核心：直接用 $e^2(n)$ 作为均方误差 $E[e^2(n)]$ 的估计值，即

$\hat{\nabla}(n) \approx \nabla[e^2(n)] = 2e(n)\nabla[e(n)] = -2e(n)X(n)$

说明: $\hat{\nabla}(n)$ 只是当个平方误差序列的梯度，而 $\nabla(n)$ 是多个平方误差序列统计平均的梯度。LMS算法就是用前者作为后者的近似，我个人理解类似用于Mini-batch版本的SGD。

最终的迭代公式为：

$W(n+1) = W(n)+2\mu e(n)X(n)$

LMS算法实践

期望信号的说明

准备两个麦克风，一个麦克风录制人+噪声，另一个麦克风远离人，只录制噪声。期望信号设定为实际的信号，而输入信号则设置为猜想的噪音信号，随着迭代的过程，输入噪音信号越接近于实际信号中包含的噪音成分。最终误差序列则为减噪后的序列。
我们已知信号为sin函数,但是不知道具体的幅值和相位。可以设置输入信号为sin函数，期望信号为实际的信号。可以通过LMS方法求的函数的幅值和相位。—陷波应用（去除单频噪声）

流程

从输入序列中 $x(n)$ 中，提取出与滤波器长度相同的一段语音段 $\hat{x}(n)$
将这段与权系数相乘
迭代误差，响应序列减去输入序列e(n)=d(n)-y(n)
根据迭代公式，更新权矢量。
⭐️将 $x(n)$ 往右平移一个数据点，重复2~5。

Matlab编程

function [yn,W,en]=LMSfilter(xn,dn,M,mu)
% 输入参数:
%     xn   输入的信号序列      (列向量)
%     dn   所期望的响应序列    (列向量)
%     M    滤波器的阶数        (标量)
%     mu   收敛因子(步长)      (标量)     要求大于0,小于xn的相关矩阵最大特征值的倒数   
% 输出参数:
%     W    滤波器的权值矩阵     (矩阵)
%          大小为M x itr,
%     en   误差序列(itr x 1)    (列向量) 
%     yn   实际输出序列         (列向量)

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
end

应用实例

Matlab自编函数测试

% 读入信号
fle='C:\Users\wangj\Documents\我的坚果云\思维导图\语音识别\声纹识别代码\宋知用voice_box\speech_signal\bluesky1.wav';
[s, fs] = v_readwav(fle);           % 读入数据文件
s = s-mean(s);                            % 消除直流分量
s = s/max(abs(s));                        % 幅值归一
N = length(s);                            % 语音长度
time = (0:N-1)/fs;                        % 设置时间刻度
noise = randn(size(s))*0.1;               % 产生随机噪声

% LMS 算法
M = 32 ;
mu = 0.001 ; 
[yn,W,en]=LMSfilter(noise,s+noise,M,mu);
subplot(3,1,1);plot(time,s+noise);xlabel('加噪后数据')
subplot(3,1,2);plot(time,en);xlabel('误差数据（原始数据）')
subplot(3,1,3);plot(time,yn);xlabel('期望数据（噪声数据）')
set(gcf,'color','w')

实验结果：

可以发现，自编的算法效果并不是特别好。下面用Matlab库里面的函数再测试一下：

Matlab自带函数测试

%%
clc;clear all;close all;

%% 产生信号源
[s,Fs] = v_readwav('C:\Users\wangj\Documents\我的坚果云\思维导图\语音识别\声纹识别代码\宋知用voice_box\speech_signal\bluesky1.wav','n');
n = length(s);
t=(0:n-1);
figure(1);
subplot(4,1,1);
plot(t,s);grid;ylabel('幅度');xlabel('时间');title('原始音频信号');
% 产生均值为0方差为0.1的噪声信号
noise = sqrt(0.1)*randn(n,1);
subplot(4,1,2);
plot(t,noise);grid;ylabel('幅度');xlabel('时间');title('相关噪声信号');
 
% 产生期望信号
dn = s + noise; % 
% audiowrite('含噪音频.wav',dn,Fs); %创建含噪音频
subplot(4,1,3);plot(t,dn);grid;
ylabel('幅度');xlabel('时间');title('含噪音频信号d(n)');

%% Matlab 自带函数
% 先设置一个训练器
lms = dsp.LMSFilter; 
[z,elms,w] = lms(noise,s+noise); 

%% 绘图部分
% 绘制去噪后的语音信号
subplot(4,1,4);plot(t,z);grid;
ylabel('幅度');xlabel('时间');title('输出序列y(n)');

%剩余噪声
figure(2);
subplot(2,1,1);
plot(t,elms);grid;ylabel('幅度');xlabel('时间');title('误差序列e(n)=d(n)-y(n)');
 
% 一小段三个信号比较
subplot(2,1,2)  ;
t= 1: length(z) ; 
plot(t,z,'r',t,elms,'g',t,s,'b') ;ylabel('幅度');xlabel('时间');
legend('去噪后的语音信号','剩余噪声','原始音频');title('一小段三个信号比较');

观察语谱图

可以发现打圈的地方，降噪的结果并不是很好。

%% 语谱图检测减噪效果：
figure(3);
wlen=200; 
window = hanning(wlen);
inc=100;
nfft = 1024 ;
subplot(3,1,1);
spectrogram(s,window,inc,nfft,Fs,'yaxis');
xlabel('纯净语音的语谱图')
subplot(3,1,2);
spectrogram(dn,window,inc,nfft,Fs,'yaxis');
xlabel('加噪语音的语谱图')
subplot(3,1,3);
spectrogram(elms,window,inc,nfft,Fs,'yaxis');
xlabel('减噪后的语谱图')

以下是LMS默认的训练超参数：

lms:
                   Method: 'LMS'
                   Length: 32
           StepSizeSource: 'Property'
                 StepSize: 0.1000
            LeakageFactor: 1
        InitialConditions: 0
           AdaptInputPort: false
    WeightsResetInputPort: false
            WeightsOutput: 'Last'

可以通过下面的代码对超参数进行设置：

lms = dsp.LMSFilter('Length',32,'Method','LMS', ...
     'StepSize',0.0001,'LeakageFactor',0.99999, ...
     'WeightsOutput','All'); % full Weights history

2. 归一化LMS滤波

通过LMS算法中对权矢量的迭代公式可以发现：

$W(n+1) = W(n)+2\mu e(n)X(n)$

$X(n)$ 的值影响到了权矢量的迭代量。当 $X(n)$ 较大时，公式可能无法收敛。

还可以发现LMS算法中的步长为一给定的常数 $\mu$ ，NLMS算法中的步长为根据时间变换的量，即

$\mu(n) = \frac{\alpha}{\|X(n)\|^2}$

考虑到分母为0等情况，最终的NLMS的迭代公式为：

$W(n+1) = W(n) + \frac{\mu}{c+\|X(n)\|_2^2}e(n)X(n)$

Matlab代码中有关dsp.LMSfilter的算法说明

无论对于不相关数据还是相关数据NLMS要比标准LMS可能呈现更快的收敛速度。

NLMS 公式推导

LMS公式存在的问题： $W(n+1)-W(n)$ 这个值很飘，受很多因素的影响。现希望自适应滤波器的权向量以最小方式改变。

原先的LMS方程为：

$e(n) = d(n)-W^T(n)X(n)$

注：这里 $W(n)$ 代表第n个时刻的权重值。而非权重系数有n个！！！

按照剧本应该求误差的均方误差最小【即误差平方的2范数最小: $min\|e^2(n)\|$ 】。

现在要求：

$d(n) = W^T(n+1)X(n)$

可以发现LMS算法的方程变为：

$e(n) = (W^T(n+1)-W^T(n))X(n)$

让 $e(n)$ 最小，等价于让 $\delta W(n+1) = W(n+1)-W(n)$ 变化最小，这里选择 $\delta W(n+1)$ 的二范数为代价函数：

$J(n) = \|\delta W(n+1)\|^2 + \lambda \cdot (d(n)-W(n+1)X(n))$

注：这里还用了拉格朗日定理。

求导等于0：

$\frac{\partial J(n)}{\partial W(n+1)} =2(W(n+1)-W(n))-\lambda\cdot X(n) = 0$

得n+1时刻的权矢量的迭代公式：

$W(n+1) = W(n) + \frac{1}{2} \lambda*X(n)$

进一步地，将上述代入限制条件，可以求解出 $\lambda$

$\begin{aligned} d(n) = W(n+1)^TX(n) = &W(n)^TX(n)+\frac{1}{2}\lambda\|X(n)\|^2 \\ = & y(n)+\frac{1}{2}\lambda\|X(n)\|^2 \end{aligned}$

进一步整理可有：

$\begin{aligned} e(n) =& d(n) - y(n) = \frac{1}{2}\lambda\|X(n)\|^2 \\ \lambda =& \frac{2e(n)}{\|X(n)\|^2} \end{aligned}$

综上，可得NLMS的权矢量的迭代公式：

$\begin{aligned} W(n+1) = &W(n) + \frac{1}{2} \frac{2e(n)}{\|X(n)\|^2}X(n) \\ =&W(n)+\frac{1}{\|X(n)\|^2}e(n)X(n) \end{aligned}$

增加一个控制增量变量的量 $\mu$ ,以及一个防止分母为0的量c（极小值），有：

$W(n+1) = W(n) + \frac{\mu}{c+\|X(n)\|_2^2}e(n)X(n)$

3. 自适应陷波

什么是陷波？

陷波器：对特定频率的信号有着很强的衰减的滤波器，也即阻带带宽极窄的带阻滤波器。

使用前提

我们需要知道原始信号里的干扰信号频率是多少时（例如最常见的50Hz工频干扰），这时我们只需要知道这个干扰信号的相位和幅度，然后就可以完全的“再现”这个干扰信号。

自适应陷波器

假设信号中的噪声是单色干扰（频率为 $\omega_0$ ）,则原始信号为： $s(t)+Acos(\omega_0 t+\phi)$ ，其中 $Acos(\omega_0 t+\phi)$ 是噪声部分。如果我们希望去除这部分的噪声，就需要知道噪声信号对应的相位和幅度。

构造输出信号 $X(n)=[x_1(n),x_2(n)]$ 为： $x_1 = cos(\omega_0 n)$ 和 $x_2 = sin(\omega_0 n)$

应用实例

宋知用案例

% pr7_1_2 
clear all; clc; close all;

%% 读入数据以及创造工频噪声
filedir=[];                             % 指定文件路径
filename='bluesky1.wav';                % 指定文件名
fle = [filedir filename];               % 构成路径和文件名的字符串
[s,fs]=wavread(fle);                    % 读入数据文件
s=s/max(abs(s));                        % 幅值归一化
N=length(s);                            % 求出信号长度
time=(0:N-1)/fs;                        % 设置时间
ns=0.5*cos(2*pi*50*time);               % 计算出50Hz工频信号
x=s+ns';                                % 语音信号和50Hz工频信号叠加

%% 设置参考输入信号
% 这部分就需要自己编写LMS算法了，不能调用Matlab自带的dsp.LMSfilter!!!

x1=cos(2*pi*50*time);                   % 设置x1和x2
x2=sin(2*pi*50*time);
w1=0.1;                                 % 初始化w1和w2
w2=0.1;
e=zeros(1, N);                          % 初始化e和y
y=zeros(1, N);
mu=0.05;                                % 设置mu
for i=1: N                              % LMS自适应陷波器滤波
  y(i)=w1 * x1(i)+ w2 * x2(i);          % 按式(7-1-29)计算y
  e(i) =x(i)-y(i);                      % 按式(7-1-30)计算e
  w1=w1+mu * e(i) * x1(i);              % 按式(7-1-31)调整w
  w2=w2+mu * e(i) * x2(i);
end
output=e';                              % 陷波器输出

%% 计算减噪前后的信噪比
snr1=SNR_singlech(s,x);                 % 计算叠加50Hz工频信号后的信噪比
snr2=SNR_singlech(s,output);            % 计算滤波后的信噪比 
snr=snr2-snr1;
fprintf('snr1=%5.4f   snr2=%5.4f    snr=%5.4f\n',snr1,snr2,snr);

%% 作图
figure(1); 
subplot 311; plot(time,s,'k'); 
ylim([-1 1 ]); title('原始语音信号');
subplot 312; plot(time,x,'k'); 
ylim([-1 1 ]); title('带噪语音信号');
subplot 313; plot(time,output,'k'); 
ylim([-1 1 ]); title('LMS滤波输出语音信号');
xlabel('时间/s')
set(gcf,'color','w')

%% 绘制频谱图
nfft = length(s) ; 
freq=linspace(0,fs/2,nfft/2+1) ;  % 设置频率刻度

% fft各个信号
fft_s = fft(s,nfft);
fft_x = fft(x,nfft);
fft_output = fft(output,nfft);

figure(2);
subplot 311; plot(freq,abs(fft_s(1:nfft/2+1))) ;
xlabel('原始语音信号');
subplot 312; plot(freq,abs(fft_x(1:nfft/2+1))) ;ylim([0,600])
xlabel('带噪语音信号');
subplot 313;  plot(freq,abs(fft_output(1:nfft/2+1))) ;
xlabel('LMS滤波输出语音信号');
set(gcf,'color','w')

%% 绘制频率图
% 设置STFT参数
wlen=200; 
window = hanning(wlen);
inc=100;

figure(3);
subplot 311; spectrogram(s,window,inc,nfft,fs,'yaxis');
xlabel('原始语音信号');
subplot 312; spectrogram(x,window,inc,nfft,fs,'yaxis');
xlabel('带噪语音信号');
subplot 313;  spectrogram(output,window,inc,nfft,fs,'yaxis');
xlabel('LMS滤波输出语音信号');
set(gcf,'color','w')

波形图和语谱图

可以发现语谱图上单波噪声还是很干净的，在频谱图上可以发现50Hz的频率噪声被减掉了。

博客案例 – LMS

clc;  
clear all;  
close all;  
%%  
%************************生成仿真信号**************************************  
Fs = 500;                                                     %设置采样频率  
t = 0:1/Fs:3;    
t = t';  
Size_t = size(t,1);  
F1 = 7;  
F2 = 13;  
F3 = 23;  
F4 = 50;  
SNR = -100;            %信干比 Unit:dB  
  
Signal = 10^(SNR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号  
Signal = Signal/max(Signal) ;
noise = 0.95*sin(2*pi*F4*t+pi/2);  
Signal_noise = Signal + noise;                               %加入50Hz工频干扰  
        
%%  
%*************************************************************************  
M = 2;                                    %定义FIR滤波器阶数  
Signal_Len = Size_t;                      %定义信号数据的个数  
niu = 0.1;                                  %算法调节步长控制因子  
y_out = zeros(Signal_Len,1);              %滤波器输出  
error_out = zeros(Signal_Len,1);          %误差输出              
w_out = zeros(Signal_Len,M);              %系数输出  
for i=1:Signal_Len  
    %数据输入  
    if i == 1           %如果是第一次进入  
        w = zeros(M,1); %初始化滤波器抽头系数  
        x = zeros(M,1); %初始化信号向量  
    end  
    
    d = Signal_noise(i);                     %输入新的期望信号  
    x = [sin(2*pi*F4*(i-1)/Fs)  
         cos(2*pi*F4*(i-1)/Fs)];             %输入新的信号矢量  
      
    %算法正体  
    y = x' * w;                              %计算滤波器输出  
    error = d - y;                           %计算误差  
    w_forward = w + niu * error * x;         %计算滤波器系数向量  
    %变量更替  
    w = w_forward;  
    %滤波结果存储  
    y_out(i) = y;  
    error_out(i) = error;  
    w_out(i,:) = w';  
end  
%%  
subplot(2,1,1);  
plot(t,Signal);  
title('原始信号');  
xlabel('时间t/s');  
subplot(2,1,2);  
plot(t,Signal_noise);  
title('加入干扰噪声的信号');  
xlabel('时间t/s');  
  
figure;  
subplot(2,1,1);  
plot(t,y_out);  
title('滤波器输出');  
xlabel('时间t/s');  
subplot(2,1,2);  
plot(t,error_out,'b');  
title('输出误差');  
xlabel('时间t/s');  
  
figure;  
plot(t(1:Signal_Len),w_out(:,1),'r',t(1:Signal_Len),w_out(:,2),'b');  
title('自适应滤波器系数');  
xlabel('时间t/s');

如果调整步长控制因子可以发现，滤波器系数需要更长的时间才会趋向稳定，最终的头部的降噪效果要打折扣。

4. RLS递归最小二乘

前言

若从代价函数的收敛方式来分，可以分为最小二乘法和梯度下降法。LMS为梯度下降法，但是通过上述实验可以发现，权重系数收敛是需要时间的，这点取决于步长。这导致了前几秒的滤波效果很差。递归最小二乘法（RLS）能有效的解决这个问题，收敛速度比LMS/NLMS收敛快一个数量级。

RLS的算法原理简述

详细理论推导可见文献👉【3】

RLS算法的关键在于用二乘方的时间平均准则取代最小均方准则，并按照时间进行迭代计算。如下：

LMS中为最小均方准则（集总平均）：

$J(k) = \|E[e(k)^2] \|$

RSL为时间平均（实数）

$\epsilon(n)=\sum_{k=0}^{n} e^{2}(k)=\min$

这里多考虑了一步，类似于移动加权平均，加入一个权重，让“新声音”的权重越大，即

$\min J_{n}(\boldsymbol{w})=\sum_{i=0}^{n} \lambda^{n-i}|e(i)|^{2}=\sum_{i=0}^{n} \lambda^{n-i}\left|y(i)-\boldsymbol{w}^{\mathrm{T}}(n) \boldsymbol{x}(i)\right|^{2}$

$0<\lambda\le1$ 为遗忘因子，这里只讨论平稳情况，取 $\lambda=1$ 。

算法流程

步骤一：初始化

$w(0)=0,P(0)=\delta^{-1}I(0<\delta\ll 1)$

其中 $\delta$ 为很小的正数，如1e-7。

当SNR较高时， $\delta$ 取小正数；当SNR较低时， $\delta$ 取大正数。

步骤二：迭代更新

计算中间量 $\pi(n)$ ，可减轻有限精读计算带来的问题

$\pi(n) = P(n-1)x(n)$
计算增益向量 $k(n)$

$\boldsymbol{k}(n)=\frac{\pi(n)}{\lambda+\boldsymbol{x}^{\mathrm{T}}(n) \pi(n)}$
计算先验估计误差

$\epsilon(n)=y(n)-\boldsymbol{w}^{\mathrm{T}}(n-1) \boldsymbol{x}(n)$
更新滤波器系数

$\boldsymbol{w}(n)=\boldsymbol{w}(n-1)+\boldsymbol{k}(n) \epsilon^{*}(n)$
更新逆相关矩阵

$\boldsymbol{P}(n)=\frac{1}{\lambda}\left[\boldsymbol{P}(n-1)-\boldsymbol{k}(n) \boldsymbol{x}^{\mathrm{H}}(n) \boldsymbol{P}(n-1)\right]$

函数部分

function [e,w]=rls(lambda,M,u,d,delta)
% recursive least squares,rls.
% Call:
% [e,w]=rls(lambda,M,u,d,delta)
%
% Input arguments:
% lambda = constant, (0,1]
% M = filter length, dim 1x1
% u = input signal, dim Nx1
% d = desired signal, dim Nx1
% delta = constant for initializaton, suggest 1e-7.
%
% Output arguments:
% e = estimation error, dim Nx1
% w = final filter coefficients, dim Mx1
% Step1:initialize
% 2017-4-4 14:34:33, Author: Gui
w=zeros(M,1);
P=eye(M)/delta;
u=u(:);
d=d(:);
% input signal length
N=length(u);
% error vector
e=d.';
% Step2: Loop, RLS
for n=M:N
    uvec=u(n:-1:n-M+1);
    e(n)=d(n)-w'*uvec;
    k=lambda^(-1)*P*uvec/(1+lambda^(-1)*uvec'*P*uvec);
    P=lambda^(-1)*P-lambda^(-1)*k*uvec'*P;
    w=w+k*conj(e(n));
end
end

应用实例

博客案例 – RLS

clc;clear all; close all;  
%%  
%************************生成仿真信号**************************************  
Fs = 500;                                                     %设置采样频率  
t = 0:1/Fs:3;    
t = t';  
Size_t = size(t,1);  
F1 = 7;  
F2 = 13;  
F3 = 23;  
F4 = 50;  
SIR = -100;            %信干比 Unit:dB  
  
Signal = 10^(SIR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号  
Signal = Signal/max(Signal) ;
noise = 0.95*sin(2*pi*F4*t+pi/2);  
Signal_noise = Signal + noise;    %加入50Hz工频干扰  
       
%%  
%*************************************************************************  
M = 2;                         %定义FIR滤波器阶数  
lamda = 0.90;                  %定义遗忘因子       
Signal_Len = Size_t;           %定义信号数据的个数  
I = eye(M);                    %生成对应的单位矩阵  
c = 0.0001;                      %小正数 保证矩阵P非奇异  
y_out = zeros(Signal_Len,1);  
Eta_out = zeros(Signal_Len,1);  
w_out = zeros(Signal_Len,M);  
for i=1:Signal_Len  
    %输入数据  
    if i == 1                 %如果是第一次进入  
        P_last = I/c;  
        w_last = zeros(M,1);    
    end  
    d = Signal_noise(i);                         %输入新的期望信号  
    x = [sin(2*pi*F4*(i-1)/Fs)  
         cos(2*pi*F4*(i-1)/Fs)];                 %输入新的信号矢量  
    %算法正体  
    K = (P_last * x)/(lamda + x'* P_last * x);   %计算增益矢量  
    y = x'* w_last;                          %计算FIR滤波器输出  
    Eta = d - y;                             %计算估计的误差  
    w = w_last + K * Eta;                    %计算滤波器系数矢量  
    P = (I - K * x')* P_last/lamda;          %计算误差相关矩阵  
    %变量更替  
    P_last = P;  
    w_last = w;  
    %滤波结果存储  
    y_out(i) = y;  
    Eta_out(i) = Eta;  
    w_out(i,:) = w';  
end  
%%  绘图部分
% 波形图
figure;  
subplot(3,1,1);  
plot(t,Signal);  
title('原始信号');  
xlabel('时间t/s');  
subplot(3,1,2);  
plot(t,Signal_noise);  
string = ['加入50Hz工频干扰的信号,信干比SIR = ',num2str(SIR),'dB'];  
title(string);  
xlabel('时间t/s');  
subplot(3,1,3);  
plot(t,Eta_out,'b',t,Signal,'r');  
legend('滤波后结果','原始信号')
title('输出误差');  
xlabel('时间t/s');  
set(gcf,'color','w')

% 滤波器的权系数
figure(2);
plot(t, w_out(:,1),'b',t,w_out(:,2),'k')
set(gcf,'color','w')
title('自适应滤波器系数');  
xlabel('时间t/s');

结果查看：

波形图对比	滤波器系数

通过结果可以发现，RSL基本没有收敛的过程，出道即收敛。毕竟是最小二乘法。

自适应滤波器减噪

0. 前言知识

什么是自适应滤波？

几种自适应滤波器的比较

常用滤波器算法

应用分类：

1. LMS滤波

LMS原理及推导

最小均方(LMS)算法

LMS算法实践

期望信号的说明

流程

Matlab编程

应用实例

Matlab自编函数测试

Matlab自带函数测试

2. 归一化LMS滤波

NLMS 公式推导

3. 自适应陷波

什么是陷波？

使用前提

自适应陷波器

应用实例

宋知用案例

博客案例 – LMS

4. RLS递归最小二乘

前言

RLS的算法原理简述

算法流程

步骤一：初始化

步骤二：迭代更新

函数部分

应用实例

博客案例 – RLS

参考网址：