clc
clear all
close all
n=500;
Sspk=zeros(1,n)';
e=zeros(1,n); % 用于存放误差
e1=zeros(1,n); % 用于存放实际回声与合成回声之间的误差
ep=zeros(1,n); % 用于存放N次运算累积的误差
w=randn(1,n)'; % 产生一个随机信号用作参考,然后用x来拟合这个信号,并用前两个点初始化滤波器第一组权向量
%算法
y=zeros(1,n);
[H,f]=freqz([-0.195 0.95],1000,16000); % 滤波器响应曲线
figure
plot(f,20*log10(abs(H)),'k')
P=1000; % LMS算法重复运算次数,用于评估P次运算产生的误差的平均水平
ee=zeros(P,n); % 用于存放平方差
erle=zeros(P,n);
ERLE=zeros(1,n);
steps = P;
h=waitbar(0,'计算进度');
for p=1:P
for i=1:n
if rand<0.5
Sspk(i)=1;
else
Sspk(i)=-1;
end
end
Secho=filter([-0.195 0.95],1,Sspk)+randn(1,n)'*sqrt(0.02);
L=2; % 滤波器阶数,考虑到两个信号之间没有延时,且这里只是用来分离出输入的x信号,而且误差不做要求,因此不需要设置太高的阶数
u=0.05; % 增益常数
wL=zeros(L,n); % 产生一个权向量矩阵
for i=(L+1):n % 计算权向量矩阵中第三组权向量到最后一组权向量相关的变量,根据x和e=x-y求下一个y,没有d(n)
X=Sspk(i-L:1:(i-1)); % 更新滤波器的参考矢量X(n)
y(i)=X'*wL(:,i); % 根据x计算i时刻输出信号
e(i)=Secho(i)-y(i);
wL(:,(i+1))=wL(:,i)+2*u*e(i)*X; % 更新i时刻滤波器的权向量
ee(p,i)=e(i)*e(i); % 更新平方差
erle(p,i)=10*log10(Secho(i)^2/e(i)^2);
end
waitbar(p/steps);
end
close(h);
eemean=sum(ee(:,:))/P;
erlemean=sum(erle)/P;
figure
subplot(2,1,1);
semilogy(eemean./max(eemean),'LineWidth',2);
subplot(2,1,2);
plot(erlemean)
figure
plot(Sspk)
hold on;
plot(y)
关键点是不要用一组数据重复计算,每次评估误差时应重新生成参考信号以及回声信号,这样求解样本的均方差才是有意义的,如果不更新数据,相当于对同一个样本重复求解N次,最终求解出来的样本均值还是自己,这样求平均值是没有意义的。
干扰误差不宜过高,否则需要非常多次的求解期望,才能让收敛曲线看上去像是一条线。
以下是上面程序的运行结果