变步长LMS算法实现

分享一个最近自己实现的变步长LMS自适应滤波的matlab程序,先来一个自适应FIR的原理图吧,方便以后复习:3
在这里插入图片描述
在这里插入图片描述
自适应噪声抵消器(ANC):
在这里插入图片描述
y 以最小均方误差逼近 n0 时,z = d - y 也以最小均方误差逼近有用信号 s 。

emmm……这里放了五种自适应滤波方法,具体方法研究大家自行找论文出处研究吧hhh就不在这里一一贴链接啦。
1:定步长LMS算法
在这里插入图片描述
2:归一化最小均方误差(NLMS)算法
在这里插入图片描述
3:基于箕舌线函数的变步长LMS算法
在这里插入图片描述
4:以及基于相关特性的改进箕舌线变步长LMS算法。
(1)变步长公式:
在这里插入图片描述
(2)更新权系数公式:
在这里插入图片描述
代码如下:

%% ===============================
% 变步长LMS算法----自适应滤波
close all;clear;clc;
M=20; %%滤波器阶次
%% 噪声信号产生器
N=1000;               %%信号长度和迭代次数控制
t=0:N-1;
snr = 1;
noise = 0.1*rand(1,N);
noise_pri = sqrt(10^(snr/10))*noise;   
noise_ref = 1*(2*rand(1,N)-1);

%% 产生正弦信号
r = 10;  % 延迟
Fs=N;              %% 采样频率
fb=24/Fs;            %% 信号频率1 
fr=28/Fs;            %% 信号频率2  干扰信号频率
s = cos(2*pi*fb*t);
p = cos(2*pi*fr*t);
n = noise_pri;

a1 = 1;
a2 = 6;
a3 = 0.5;

b1 = 0;
b2 = 6;
b3 = 0.5;
input_ref = b1*s + b2*p + b3*n ; %% 滤波器输入信号
d_pri = a1*s + a2*p + a3*n; %理想输出
% d_pri = input_ref(r:N);  % 可以采用延迟做输入
x = input_ref;
figure(1);
subplot(3,1,1);plot(s);grid on;title(' 有用信号 s(n) ');
subplot(3,1,2);plot(d_pri);grid on;title(' 主输入 d(n) ');
subplot(3,1,3);plot(input_ref);grid on;title(' 参考输入 x(n)');

%% 自适应滤波
g = 100;
tz = (eig(input_ref'*input_ref))';
rho_max = max(eig(input_ref'*input_ref));   % 输入信号相关矩阵的最大特征值
mu_max = 2/rho_max; % 收敛范围极限值
for type = 1:5
    for q = 1:g
%         noise = 0.1*rand(1,N); % 每次加入不同的噪声
%         noise_pri = sqrt(10^(snr/10))*noise;   
%         input_ref = b1*s + b2*p + b3*noise_pri;
%         d_pri = a1*s + a2*p + a3*noise_pri;
        
        xx=input_ref;          % 未知系统输入等于输出           
        W=zeros(1,M);     % 滤波器权向量
        W1=zeros(1,M);
        W2=zeros(1,M);
        e=zeros(1,N-r);        % 误差
        y=zeros(1,N-r);
        mu=zeros(1,N-r);
        pp = zeros(g,N-M-r);
%         bi = zeros(g,N-M-r);
        for j=M:N-r
            y(j)=W*(xx(j:-1:j-M+1))';
            e(j)=d_pri(j)-y(j);
            err(j) = (e(j) - s(j));
            % ---------------- 关于mu -------------------
            if type == 1   
                % 方法1 ----- LMS 定步长
                mu(j)=mu_max/4;    
                W = W +2*mu(j)*e(j)*(xx(j:-1:j-M+1));
            elseif type == 2
                % 方法2 ----- NLMS
%                 mu(j) = mu_max./(1+abs(e(j).^2)-0.2); 
                mu(j) = 0.015/((xx(j:-1:j-M+1))*(xx(j:-1:j-M+1))');
                W = W +2*mu(j)*e(j)*(xx(j:-1:j-M+1));
            elseif type == 3
                % 方法3 ----- 箕舌线
                b2 = 0.95*mu_max;a2 = 1;
                mu(j) = b2*(1-(1/(a2*e(j)*e(j)+1+0.000001))); 
                W = W +2*mu(j)*e(j)*(xx(j:-1:j-M+1));    
            elseif type == 4
                % 方法4 ----- 基于相关特性的改进箕舌线变步长LMS1
                b2 = 0.95*mu_max;a2 = 1;
                mu(j) = b2*(1-(1/(a2*e(j)*e(j-1)+1+0.000001))); 
                W = W +2*mu(j)*e(j)*(xx(j:-1:j-M+1));
            else
                % 方法5 --基于相关特性的改进箕舌线变步长LMS2
                b2 = 0.95*mu_max;a2 = 1;
                mu(j) = b2*(1-(1/(a2*e(j)*e(j-1)+1+0.000001))); 
                d = 0.05000;               
                W = W2 +2*mu(j)*e(j)*(xx(j:-1:j-M+1))+d*(W2-W1);
                mm = d*(W2-W1);
                W1 = W2;
                W2 = W;
            end  
%             sum = sum + err(j);
        end                   
        pp(q,:) = (err((M+1):(N-r))).^2;
    end
    bi = mean(pp,1);
%% ================================ 绘图 ===================================
    figure(2); % 绘制输出信号
    if type == 1
        subplot(521);plot(y);title(' 滤波器输出 y(LMS)');grid on;
        yn_fft = abs(fft(y));
        subplot(522);stem(yn_fft(1:50));title(' 滤波器输出 y 频谱');grid on;
    elseif type == 2 
        subplot(523);plot(y);title(' 滤波器输出 y(NLMS)');grid on;
        yn_fft = abs(fft(y));
        subplot(524);stem(yn_fft(1:50));title(' 滤波器输出 y 频谱');grid on;
    elseif type == 3
        subplot(525);plot(y);title(' 滤波器输出 y(箕舌线)');grid on;
        yn_fft = abs(fft(y));
        subplot(526);stem(yn_fft(1:50));title(' 滤波器输出 y 频谱');grid on;
    elseif type == 4 
        subplot(527);plot(y);title(' 滤波器输出 y(改进1)');grid on;
        yn_fft = abs(fft(y));
        subplot(528);stem(yn_fft(1:50));title(' 滤波器输出 y 频谱');grid on;
    else 
        subplot(529);plot(y);title(' 滤波器输出 y(改进2)');grid on;
        yn_fft = abs(fft(y));
        subplot(5,2,10);stem(yn_fft(1:50));title(' 滤波器输出 y 频谱');grid on;
    end
    e_fft = abs(fft(e));
    figure(3);% 绘制误差信号
    if type == 1
        subplot(521);plot(e);title(' 输出 e (LMS)');grid on;
        subplot(522);plot(e_fft);title(' 频谱 (LMS)');grid on;
    elseif type == 2 
        subplot(523);plot(e);title(' 输出 e (NLMS)');grid on;
        subplot(524);plot(e_fft);title(' 频谱 (NLMS)');grid on;
    elseif type == 3
        subplot(525);plot(e);title(' 输出 e (箕舌线)');grid on;
        subplot(526);plot(e_fft);title(' 频谱 (箕舌线)');grid on;
    elseif type == 4 
        subplot(527);plot(e);title(' 输出 e (改进1)');grid on;
        subplot(528);plot(e_fft);title('  频谱 (改进1)');grid on;
    else 
        subplot(529);plot(e);title(' 输出 e (改进2)');grid on;
        subplot(5,2,10);plot(e_fft);title(' 频谱 (改进2)');grid on;
    end 
    
     figure(4);% 绘制误差信号
    if type == 1
        subplot(521);plot(err);title(' 误差信号 err (LMS)');grid on;
        subplot(522);plot(bi);title(' 均方误差 (LMS)');grid on;
    elseif type == 2 
        subplot(523);plot(err);title(' 误差信号 err (NLMS)');grid on;
        subplot(524);plot(bi);title(' 均方误差 (NLMS)');grid on;
    elseif type == 3
        bi3 = bi;
        subplot(525);plot(err);title(' 误差信号 err(箕舌线)');grid on;
        subplot(526);plot(bi);title(' 均方误差 (箕舌线)');grid on;
    elseif type == 4 
        subplot(527);plot(err);title(' 误差信号 err (改进1)');grid on;
        subplot(528);plot(bi);title(' 均方误差 (改进1)');grid on;
        bi4 = bi;
    else
        bi5 = bi;
        subplot(529);plot(err);title(' 误差信号 err (改进2)');grid on;
        subplot(5,2,10);plot(bi);title(' 均方误差 (改进2)');grid on;
        
    end 
%     figure(5)
%     plot((mu));hold on;legend('LMS','NLMS','箕舌线','改进1','改进2');
end
% "学习"曲线是表示均方误差对时间或迭代次数的关系曲线.
% 均方误差是权系数矢量的二次型.由于权系数的过渡过程是指数函数,
% 指数函数的平方仍然为指数函数,但其时间常数是原来的一半,所以"学习"曲线应具有指数特性

实验效果图 :
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
有问题欢迎大家批评指正:3

  • 15
    点赞
  • 123
    收藏
  • 打赏
    打赏
  • 12
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:大白 设计师:CSDN官方博客 返回首页
评论 12

打赏作者

Sakura.66

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值