function en=LMS(xn,dn,M,mu,itr)
% 输入参数:
% xn 输入的信号序列 (列向量)
% dn 所期望的响应序列 (列向量)
% M 滤波器阶数 (标量)
% mu 步长因子(0~1) (标量)
% itr 迭代长度 (标量) 默认为xn的长度,M<itr<length(xn)
%中间参数:
% W 滤波器的权值矩阵 (矩阵)
% 大小为M 行× itr 列
% yn 实际输出序列 (列向量)
% 输出参数:
% en 误差序列(itr x 1) (列向量)
% 参数个数必须为4个或5个
%4个时,即以xn长度为迭代次数;5个时,即以给出的第五个参数为迭代次数
if nargin == 4 % 4个时递归迭代的次数为xn的长度
itr = length(xn);
elseif nargin == 5 % 5个时满足M<itr<length(xn)
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处,迭代才能开始)
x=xn(k:-1:k-M+1); % 滤波器M个抽头的输入(与论文同,反向回溯)
y=W(:,k-1).'*x; % 滤波器的输出(权向量(列向量)之转置乘输入(列向量),得到一个标量)
en(k)=dn(k)-y ; % 第k次迭代的误差
tr=2*sum(eig(x*x.'));
W(:,k) =W(:,k-1)+2*mu*en(k)*x/tr;
end
% 求最优时滤波器的输出序列
yn = inf * ones(size(xn)); %产生长度为xn大小的列向量,初始填充为无穷(对应实际的物理意义是,滤波器尚未开始工作,次级声源尚未开启)
for k=M:itr
x=xn(k:-1:k-M+1);
yn(k)=W(:,k-1).'* x;
end
for k=itr+1:length(xn)
x=xn(k:-1:k-M+1);
yn(k)=W(:,end).'* x;
end
en=dn-yn;
end