自适应滤波器和主动降噪——LMS、NLMS、RLS
- LMS滤波算法
% ------------------------------------------------------------------------
% myLMS.m - least mean squares algorithm
%
% Usage: [e, y, w] = myLMS(d, x, mu, M)
%
% Inputs:
% d - the vector of desired signal samples of size Ns, 参考信号
% x - the vector of input signal samples of size Ns, 输入信号
% mu - the stepsize parameter, 步长
% M - the number of taps. 滤波器阶数
%
% Outputs:
% e - the output error vector of size Ns
% y = output coefficients
% w - filter parameters
%
% ------------------------------------------------------------------------
function [e, y, w] = myLMS(d, x, mu, M)
Ns = length(d);
if (Ns <= M)
print('error: 信号长度小于滤波器阶数！');
return;
end
if (Ns ~= length(x))
print('error: 输入信号和参考信号长度不同！');
return;
end
x = x;
xx = zeros(M,1);
w1 = zeros(M,1);
y = zeros(Ns,1);
e = zeros(Ns,1);
for n = 1:Ns
xx = [xx(2:M);x(n)];
y(n) = w1' * xx;
e(n) = d(n) - y(n);
w1 = w1 + mu * e(n) * xx;
w(:,n) = w1;
end
end
- NLMS滤波算法
% ------------------------------------------------------------------------
% myNLMS.m - Normalized least mean squares algorithm
%
% Usage: [e, y, w] = myLMS(d, x, mu, M, a)
%
% Inputs:
% d - the vector of desired signal samples of size Ns, 参考信号
% x - the vector of input signal samples of size Ns, 输入信号
% mu - the stepsize parameter, 步长
% a - the bias parameter, 偏置参数
% M - the number of taps. 滤波器阶数
%
% Outputs:
% e - the output error vector of size Ns
% y - output coefficients
% w - filter parameters
%
% ------------------------------------------------------------------------
function [e,y, w] = myNLMS(d, x, mu, M, a)
Ns = length(d);
if (Ns <= M)
print('error: 信号长度小于滤波器阶数！');
return;
end
if (Ns ~= length(x))
print('error: 输入信号和参考信号长度不同！');
return;
end
x = x;
xx = zeros( M,1);
w1 = zeros( M,1);
y = zeros(Ns,1);
e = zeros(Ns,1);
for n = 1:Ns
xx = [xx(2:M);x(n)];
y(n) = w1' * xx;
k = mu/(a + xx'*xx);
e(n) = d(n) - y(n);
w1 = w1 + k * e(n) * xx;
w(:,n) = w1;
end
end
- RLS滤波算法
% ------------------------------------------------------------------------
% myRLS.m - least mean squares algorithm
%
% Usage: [e, y, w] = myRLS(d, x, lamda, M)
%
% Inputs:
% d - the vector of desired signal samples of size Ns, 参考信号
% x - the vector of input signal samples of size Ns, 输入信号
% lamda - the weight parameter, 权重
% M - the number of taps. 滤波器阶数
%
% Outputs:
% e - the output error vector of size Ns
% y - output coefficients
% w - filter parameters
%
% ------------------------------------------------------------------------
function [e, y, w] = myRLS(d, x,lamda,M)
Ns = length(d);
if (Ns <= M)
print('error: 信号长度小于滤波器阶数！');
return;
end
if (Ns ~= length(x))
print('error: 输入信号和参考信号长度不同！');
return;
end
I = eye(M);
a = 0.01;
p = a * I;
x = x; %在输入信号x前补上M-1个0，使输出y与输入具有相同长度
w1 = zeros(M,1);
y = zeros(Ns, 1);
e = zeros(Ns, 1);
xx = zeros(M,1);
for n = 1:Ns
xx = [x(n); xx(1:M-1)];
k = (p * xx) ./ (lamda + xx' * p * xx);
y(n) = xx'*w1;
e(n) = d(n) - y(n);
w1 = w1 + k * e(n);
p = (p - k * xx' * p) ./ lamda;
w(:,n) = w1;
end
end
- 自适应滤波器demo.
% ------------------------------------------------------------------------
%
% demo1.m - adaptive filter demo（自适应滤波算法应用demo）
% Including LMS、NLMS、RLS algorithm
% Including：
% 1、echo cancellation 音频回声消除
% 2、audio + white noise 音频白噪声消除
% 3、audio + single frequency noise 音频+单频噪声消除
% 4、single frequency signal + white noise 单频信号+白噪声消除
% 5、multi-frequency signal + single frequency signal noise 多频信号+单频声消除
% Parameters：
% x ： input signal 输入信号
% d ： reference signal 参考信号
% y ： output signal 输出信号
% e ： error signal 误差信号
% mu ： LMS stepsize LMS算法步长
% mu2 ： NLMS stepsize NLMS算法步长
% a ： NLMS bias NLMS算法偏置参数
% lamda ： RLS weight RLS算法权重
%
% ------------------------------------------------------------------------
close all;clear;clc;
%% 1、echo cancellation (音频回声)
% [d, fs] = audioread('handel.wav');
% x = audioread('handel_echo.wav');
% mu = 0.05;
% mu2 = 0.05;
% a = 0.1;
% lamda = 0.99;
% M = 40;
%% 2、audio + white noise (音频+白噪声)
% [d,fs] = audioread('handel.wav');
% n = length(d);
% noise = wgn(1, n, -20)';
% x = d + noise;
% mu = 0.001;
% mu2 = 0.1;
% a = 0.01;
% lamda = 0.99;
% M = 80;
%% 3、audio + single frequency noise (音频+单频噪声)
% [d, fs] = audioread('handel.wav');
% n = length(d);
% T = n/fs;
% t = 0:1/fs:T-1/fs;
% noise = cos(2*pi*t*370)';
% x = d + noise;
%
% mu = 0.005;
% mu2 = 0.05;
% a = 0.1;
% lamda = 0.9999;
% M = 40;
%% 4、single frequency signal + white noise (单频+白噪声)
fs = 8000;
t = 0:1/fs:5;
noise = wgn(1, length(t),-20)';
d = cos(2*pi*t*261.625)';
x = noise + cos(2*pi*t*261.625)';
mu = 0.0001;
mu2 = 0.001;
a = 0.1;
lamda = 0.999;
M = 40;
%% 5、multi-frequency signal + single frequency signal noise (多频+单频噪声)
% fs = 16000;
% t = 0:1/fs:4;
% d = cos(2*pi*t*200)';
% x = cos(2*pi*t*100)' + d;
% mu = 0.01;
% mu2 = 0.1;
% a = 0.1;
% lamda = 0.9999;
% M = 50;
%% LMS\NLMS\RLS performance (LMS\NLMS\RLS性能比较)
% run algorithm (运行算法并计算时间)
tic
[e1, y1, w1] = myLMS(d, x, mu, M);
toc
tic
[e2, y2, w2] = myNLMS(d, x,mu2, M, a);
toc
tic
[e3, y3, w3] = myRLS(d, x,lamda,M);
toc
% 画出输入信号、参考信号、滤波输出、误差
figure()
subplot(4,2,1)
plot([1:length(x)]/fs,x);
xlabel('time');
title('x(n)');
subplot(4,2,2)
plot([1:length(d)]/fs,d);
xlabel('time');
title('d(n)');
subplot(4,2,3)
plot([1:length(y1)]/fs,y1);
xlabel('time');
title('LMS y(n)');
subplot(4,2,5)
plot([1:length(y2)]/fs,y2);
xlabel('time');
title('NLMS y(n)');
subplot(4,2,7)
plot([1:length(y3)]/fs,y3);
xlabel('time');
title('RLS y(n)');
subplot(4,2,4)
plot([1:length(e1)]/fs,e1);
xlabel('time');
title('LMS e(n)');
subplot(4,2,6)
plot([1:length(e2)]/fs,e2);
xlabel('time');
title('NLMS e(n)');
subplot(4,2,8)
plot([1:length(e3)]/fs,e3);
xlabel('time');
title('RLS e(n)');
- 主动噪声控制demo.
% ------------------------------------------------------------------------
%
% demo2.m - ANC demo（使用自适应滤波算法的ANCdemo）
% Including LMS、NLMS、RLS algorithm
% Including：
% 1、audio + white noise 音频白噪声消除
% 2、single frequency signal + white noise 单频白噪声消除
% Parameters：
% x ： input signal 输入信号
% d ： reference signal 参考信号
% y ： output signal 输出信号
% e ： error signal 误差信号
% mu ： LMS stepsize LMS算法步长
% mu2 ： NLMS stepsize NLMS算法步长
% a ： NLMS bias NLMS算法偏置参数
% lamda ： RLS weight RLS算法权重
%
% System：
% signal+noise_____________d(n)___________
% +↓
% noise’—x(n)—【filter】——y(n)— - —-O——e(n)——
% ↑_____________________|
%
% ------------------------------------------------------------------------
close all;clear;clc;
%% 1、audio + white noise（音频+白噪声）
[signal,fs] = audioread('handel.wav');
noise = wgn(length(signal), 1, -20);
d = signal + noise;
x = sin(1./(1+exp(-noise)));
mu = 0.1;
mu2 = 0.5;
a = 0.01;
lamda = 0.999;
M = 20;
%% 2、single frequency signal + white noise（单频+白噪声）
% fs = 8000;
% t = 0:1/fs:4;
% signal = cos(2*pi*t*20)';
% noise = wgn(1,length(t),-20)';
% d = noise + signal;
% x = sin(1./(1+exp(-noise)));
%
% mu = 0.1;
% mu2 = 0.8;
% a = 0.01;
% lamda = 0.9999;
% M = 20;
%% LMS\NLMS\RLS performance（LMS\NLMS\RLS性能比较）
% run algorithm (运行算法)
tic
[e1, y1, w1] = myLMS(d, x, mu, M);
toc
tic
[e2, y2, w2] = myNLMS(d, x,mu2, M, a);
toc
tic
[e3, y3, w3] = myRLS(d, x,lamda,M);
toc
% 画出输入信号、参考信号、滤波输出、误差
figure()
subplot(4,2,1)
plot([1:length(x)]/fs,x);
xlabel('time');
title('x(n)');
subplot(4,2,2)
plot([1:length(d)]/fs,d);
xlabel('time');
title('d(n)');
subplot(4,2,3)
plot([1:length(y1)]/fs,y1);
xlabel('time');
title('LMS y(n)');
subplot(4,2,5)
plot([1:length(y2)]/fs,y2);
xlabel('time');
title('NLMS y(n)');
subplot(4,2,7)
plot([1:length(y3)]/fs,y3);
xlabel('time');
title('RLS y(n)');
subplot(4,2,4)
plot([1:length(e1)]/fs,e1);
xlabel('time');
title('LMS e(n)');
subplot(4,2,6)
plot([1:length(e2)]/fs,e2);
xlabel('time');
title('NLMS e(n)');
subplot(4,2,8)
plot([1:length(e3)]/fs,e3);
xlabel('time');
title('RLS e(n)');
% 画出参考信号与滤波输出的差值(ANC输出信号)
figure()
subplot(3,1,1)
plot(e1)
title('LMS ANC输出')
subplot(3,1,2)
plot(e2)
title('NLMS ANC输出')
subplot(3,1,3)
plot(e3)
title('RLS ANC输出')
% 比较稳定后的信噪比
% xx1 = clearspeech(length(x)-3000:length(x));
% ee1 = e1(length(x)-3000:length(x));
% ee2 = e2(length(x)-3000:length(x));
% ee3 = e3(length(x)-3000:length(x));
% yy1 = y1(length(x)-3000:length(x));
% yy2 = y2(length(x)-3000:length(x));
% yy3 = y3(length(x)-3000:length(x));
% SNR1 = snr(xx1,ee1)
% SNR2 = snr(xx1,ee2)
% SNR3 = snr(xx1,ee3)
%% 试听RLS ANC的输出结果（RLS效果最好）
% sound(x-y3,fs);
- 算法性能比较demo.
% ------------------------------------------------------------------------
% demo3.m - 自适应滤波算法性能比较
% Including LMS、NLMS、RLS algorithm
% Including：
% 1、不同步长的滤波器参数更新曲线对比
% 2、相同步长的LMS算法和NLMS算法滤波器权重更新曲线对比
% Parameters：
% x ： input signal 输入信号
% d ： reference signal 参考信号
% y ： output signal 输出信号
% e ： error signal 误差信号
% mu ： LMS stepsize LMS算法步长
% mu2 ： NLMS stepsize NLMS算法步长
% a ： NLMS bias NLMS算法偏置参数
% lamda ： RLS weight RLS算法权重
%
% ------------------------------------------------------------------------
close all;clear;clc;
%% audio + single frequency noise （音频+单频噪声）
% [d,fs] = audioread('handel.wav');
% n = length(d);
% T = n/fs;
% t = 0:1/fs:T-1/fs;
% noise = cos(2*pi*t*700)';
% x = d + noise;
%% single frequency signal + white noise （单频+白噪声）
fs = 8000;
t = 0:1/fs:2;
noise = wgn(length(t),1,-20);
d = cos(2*pi*t*723)' + sin(2*pi*t*456)';
x = noise + d;
%% LMS\NLMS\RLS性能比较
% set parameters (设置参数)
mu = [0.0005 0.001 0.005 0.01];
mu2 = [0.005 0.01 0.05 0.1];
a = 0.01;
lamda = [1 0.9999 0.999 0.99];
M = 20;
% 画出滤波器最后一个参数的变化曲线
figure(1);
for i = 1:length(mu)
[e1, ~, w1] = myLMS(d, x, mu(i), M);
c1(i) = {['\mu = ',num2str(mu(i))]};
subplot(2,1,1)
plot(w1(M,:)','LineWidth', 1)
hold on
subplot(2,1,2)
plot(e1)
hold on
end
subplot(2,1,1)
legend(c1)
title('不同步长的LMS算法滤波器参数收敛情况,M=20')
subplot(2,1,2)
legend(c1)
title('不同步长的LMS算法误差收敛情况,M=20')
figure(2)
for i = 1:length(mu2)
[e2, ~, w2] = myNLMS(d, x,mu2(i), M, a);
c2(i) = {['\mu = ',num2str(mu2(i))]};
subplot(2,1,1)
plot(w2(M,:)','LineWidth', 1)
hold on
subplot(2,1,2)
plot(e2)
hold on
end
subplot(2,1,1)
legend(c2)
title('不同步长的NLMS算法滤波器参数收敛情况,M=20')
subplot(2,1,2)
legend(c2)
title('不同步长的NLMS算法误差收敛情况,M=20')
figure(3)
for i = 1:length(lamda)
[e3, ~, w3] = myRLS(d, x,lamda(i),M);
c3(i) = {['\lambda = ',num2str(lamda(i))]};
subplot(2,1,1)
plot(w3(M,:)','LineWidth', 1)
hold on
subplot(2,1,2)
plot(e3)
hold on
end
subplot(2,1,1)
legend(c3)
title('不同权重的RLS算法滤波器参数收敛情况,M=20')
subplot(2,1,2)
legend(c3)
title('不同权重的RLS算法误差收敛情况,M=20')
figure(4)
[e1, ~, w1] = myLMS(d, x, 0.001, M);
[e2, ~, w2] = myNLMS(d, x, 0.01, M, a);
[e3, ~, w3] = myRLS(d, x, 0.9999,M);
subplot(2,1,1)
plot(w2(M,:)','LineWidth', 1);
hold on
plot(w1(M,:)','LineWidth', 1);
plot(w3(M,:)','LineWidth', 1);
legend({'NLMS','LMS','RLS'})
title('滤波器参数收敛情况对比：\mu1 = 0.001, \mu2 = 0.01, \lambda = 0.9999,M=20')
subplot(2,1,2)
plot(e2);
hold on
plot(e1);
plot(e3);
legend({'NLMS','LMS','RLS'})
title('误差收敛情况对比：\mu1 = 0.001, \mu2 = 0.01,\lambda = 0.9999,M=20')
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步