可以与之前的dsp.LMSFilter作为对比,从原理来实现自适应滤波器中的最小均方误差算法(Least Mean Square)
注意:本程序仅兼容2018a以后版本
这个程序有一个地方目前还存在一些问题,就是求取自相关函数R时所用到的方法,我通过尝试自行编写自相关函数求取方程和matlab自带的corrmtx函数作为对比,发现数值有所不同,应该是matlab自带的corrmtx函数和书本上的求解方式并不相同,所以使用时需要自行判断使用哪种方式。
clc
clear
%%构建信号%%
N=256;%信号长度
for j=1:50 %进行50次lms运算
t=1:N;
d=sin(0.1*pi*t);
figure(1)
plot(t,d);
axis([0,N,-2,2]);
%%加噪声%%
x=awgn(d,20);%awgn 加性高斯白噪声 信噪比20DB
% 一种简单的计算自相关矩阵方法
% R=x'*x;%计算自相关矩阵
%%直接计算自相关矩阵 设自相关矩阵长度为L自己写
% L=256;%自相关长度
% Rxx=zeros(L,L);
% for m=1:L
% for n=1:L
% mm=abs(m-n);
% temp=0;
% for nn=1:N-mm
% temp=temp+x(nn)*x(nn+mm);
% end
% Rxx(m,n)=temp/(N-mm);
% end
% end
% [V,D]=eig(Rxx);%特征值分解求u
% ma=max(max(D));
% u=1/ma;
% 一维实值信号x的自相关矩阵Rxx应为实对称的toeplitz矩阵,
%而一维实值信号x,y的互相关矩阵Rxy为非对称的toeplitz阵,
%matlab提供的corrmtx函数产生的并非通常意义下的autocorrelation matrix
m= 255;% dfine the time lag m+1,and m+1<=n
[rx,r]=corrmtx(x,m);%corrmtx为自相关矩阵
R=rx'*rx;
[V,D]=eig(R);%提取自相关矩阵的特征值对角阵D
maxlamuda=max(D,[],'all');%提取自相关矩阵最大特征值
u=1/maxlamuda-0.01;%保证LMS算法收敛并保持稳定的条件
%%LMS迭代滤波%%
T=50;
k=25;
p=zeros(T,N-k);%50次 误差矩阵
e=zeros(1,N);%误差矩阵
w=zeros(1,k);%滤波器长度系数
y=zeros(1,N);%滤波输出
y(1:k)=x(1:k);%假设滤波器前25点为输入值
for i=(k+1):N%开始迭代更新参数
Xn=x(i:-1:(i-k+1));%输入
y(i)=w*Xn';%输出
e(i)=d(i)-y(i);%误差
w=w+2*u*e(i)*Xn;%权值更新
end
p(j,:)=(e(k+1:N)).^2;%误差信号功率
end
for i=1:N-k
E(i)=sum(p(:,i))/T;%50次误差取平均
end
figure(2)
t=1:N-k;
plot(t,E);
Copyright © 2020 by RichardYang. All rights reserved.
仅供参考,严禁转载,感谢。
算法工作模式
参数:
- 抽头权系数向量w(n)
- 输入向量x(n)
- 期望输出d(n)
- 滤波器输出y(n)
- 滤波器抽头权更新参数w(n+1)
流程:
- 滤波: y ( n ) = ω T ( n ) x ( n ) \displaystyle y(n)={{\omega }^{T}}(n)x(n) y(n)=ωT(n)x(n)
- 误差估计: e ( n ) = d ( n ) − y ( n ) \displaystyle e(n)=d(n)-y(n) e(n)=d(n)−y(n)
- 抽头权向量更新: ω ( n + 1 ) = ω ( n ) + 2 μ e ( n ) x ( n ) \displaystyle \omega (n+1)=\omega (n)+2\mu e(n)x(n) ω(n+1)=ω(n)+2μe(n)x(n)
注意:
- 为了保证算法的收敛和稳定性必须有: 0 < μ < 1 λ max \displaystyle 0<\mu <\frac{1}{{{{\lambda }_{{\max }}}}} 0<μ<λmax1
- λ max \displaystyle {{\lambda }_{{\max }}} λmax为输入自相关矩阵R的最大特征值(公式证明,但实际使用时可能会出现问题,通常不这么求最大步长,但测试可以使用 1)
自适应滤波器原理及Matlab仿真应用(原书第2版) ↩︎