湖南大学计算机与通信学院
课程作业2
题 目:基于LMS和RLS的自适应滤波器的应
用仿真
基于LMS和RLS的自适应滤波器应用仿真
1. 自适应滤波原理
自适应滤波器是指利用前一时刻的结果,自动调节当前时刻的滤波器参数,以适应信号和噪声未知或随机变化的特性,得到有效的输出,主要由参数可调的 数字滤波器和自适应算法两部分组成,如图1所示
图1 自适应滤波器原理图
x(n)称为输入信号,y(n)称为输出信号,d(n)称为期望信号或者训练信号,e(n)为误差僖号,其中,e(n)=d(n)-y(n).自适应滤波器的系数(权值)根据误差信号e(n),通过一定的自适应算法不断的进行改变,以达到使输出信号y(n)最接近期望信号
图中参数可调的数字滤波器和自适应算法组成自适应滤波器。自适应滤波算法是滤波器系数权值更新的控制算法,根据输入信号与期望信号以及它们之间的误差信号,自适应滤波算法依据算法准则对滤波器的系数权值进行更新,使其能够使滤波器的输出趋向于期望信号。 原理
记数字滤波器脉冲响应为:
T
h(k)=[h0(k) h1(k) … hn-1(k)]输入采样信号为:
x(k)=[x(k) x(k-1) … x(k-n-1)] 误差信号为:
e(k)?y(k)?y(k) e(k)?y(k)?h(k)x(k)
T^优化过程就是最小化性能指标J(k),它是误差的平方和:
J(k)??[y(i)?hT(k)x(i)]2
i?1k求使J(k)最小的系数向量h(k),即使J(k)对h(k)的导数为零,也就是J(k)的表达式代入,得:
dJ(k)?0。把dh(k)
2?[y(i)?hi?1TkT(k)x(i)]x(i)?0
和
?xi?1k(i)y(i)?h(k)?x(i)xT(i)
Ti?1k由此得出滤波器系数的最优向量:
h(k)?T?xi?1ki?1kT(i)y(i)
T?x(i)x(i)这个表达式由输入信号自相关矩阵cxx(x)和输入信号与参考信号的相关矩阵
cyx(k)组成,如下所示,维数都为(n,n):
cxx(k)??xi?1kkT(i)x(i)
cyx(k)??xi?1T(i)y(i)
系数最优向量也可以写成如下形式:
?1hT(k)opt?cyx(k)cxx(k)
自相关和互相关矩阵的递归表达式如下:
cxx(k)?cxx(k?1)?x(k)xT(k) cyx(k)?cyx(k?1)?y(k)x(k)
T把cyx(k)的递归表达式代入系数向量表达式,得:
?1 hT(k)?cyx(k)cxx(k)
即
h(k)?[cyx(k?1)?x(k)y(k)]cxx(k)
TT?1考虑到
cyx(k?1)?h(k?1)cxx(k?1)
T可以记
h(k)?cxx(x)[cxx(k?1)h(k?1)?y(k)x(k)]
?1
用前面得到的表达式求出cxx(k?1),并代入上式:
?1 h(k)?cxx(x){[cxx(k)?x(k)xT(k)]h(k?1)?y(k)x(k)} ?1或 h(k)?h(k?1)?cxx(x)[y(k)x(k)?x(k)xT(k)h(k?1)]
则滤波器系数的递归关系式可以记作
?1 h(k)?h(k?1)?cxx(x)[y(k)x(k)?x(k)xT(k)h(k?1)]
其中
e(k)?y(k)?xT(k)h(k?1)
e(k)表示先验误差。只因为它是由前一个采样时刻的系数算出的,在实际中,很多时候由于h(k)计算的复杂度而不能应用于实时控制。用δ,I代换cxx(k),其中:δ为自适应梯度,I为辨识矩阵(n,n) 这时
h(k)?h(k?1)??x(k)e(k)
这时就是一个最小均方准则问题。
2. LMS自适应滤波器举例
自回归过程的自适应预估器
自回归过程是用来描述伴随一些可能性规律出现的统计现象的瞬时估计的随机过程。一阶自回归模型的公式如下:
y(k)??a1y(k?1)?b(k)
a1是模型的唯一参数,b(k)是零均值白噪声。用一个自适应滤波器生成一个可以对参数a1进
行一步预测的一阶自适应预估器。LMS算法可由如下方程表示:
e(k)?y(k)?^^a1^y(k?1)
1
a(k?1)?1a(k)???y(k?1)e(k)
取N个点估计参数a1,为获取平均值重复M次。而且分别对δ=0.01,δ=0.05,δ=0.1进行计算。参数a1固定在-0.6。 程序清单如下: N=500;
M=20; n=1;
a1=-0.8;
h=zeros(M,n+1,3); e=zeros(M,n,3); for d=1:3
if d==1 delta=0.01; else delta=0.05*(d-1); end; for k=1:M
b=0.2*randn(1,N); y(1)=1; for i=2:N
y(i)=-a1*y(i-1)+b(i); end for i=n+1:N
e(k,i,d)=y(i)-h(k,i,d)*y(i-1);
h(k,i+1,d)=h(k,i,d)+delta*y(i-1)*e(k,i,d); end
end
end
for d=1:3 for i=1:N
em(i,d)=0; hm(i,d)=0; for j=1:M
em(i,d)=em(i,d)+e(j,i,d)^2; hm(i,d)=hm(i,d)+h(j,i,d); end end end
figure(1)
semilogy(1:150,em(1:150,1)),hold on semilogy(1:150,em(1:150,2),'r'),hold on semilogy(1:150,em(1:150,3),'g'),hold off axis([0 150 0.01 1]),grid title('Mean square error ') xlabel('Samples')
gtext('\\leftarrowd=0.01'); gtext('\\leftarrowd=0.05'); gtext('\\leftarrowd=0.1');
figure(2),plot(1:N,hm(1:N,1)),hold on plot(1:N,hm(1:N,2),'r'),hold on
plot(1:N,hm(1:N,3),'g'),hold off,grid title('Filter coeffcient evalution') xlabel('Samples'),