RLS算法辨识系统.doc

自适应滤波器(分别采用LMS和RLS算法)对IIR系统辨识
1. LMS和RLS算法简介
LMS自适应滤波算法是以均方误差最小为准则的，其意义是是滤波器的输出信号和需要信号的误差平方的统计平均值最小，这个准则是依据输入数据的长期统计特性寻求最佳滤波的。而最小二乘（RLS）算法是只依据一组数据来寻求最佳滤波的方法。因此，LMS滤波是统计意义上的最佳滤波器，而RLS滤波是针对不同的数据组生成不同的最佳滤波器。RLS算法实际上是FIR维纳滤波器的一种时间递归算法，它是严格以最小二乘方准则为依据的算法。它的主要优点是收敛速度快，因此，首先在快速信道均衡，实时系统辨识和时间序列分析中得到广泛应用其主要缺点是每次迭代计算量很大(对于阶横向滤波器，计算量数量级为），因此，在信号处理中它的应用曾一度收到限制。但是近年来人们重新对它产生了兴趣，主要是因为它具有收敛速度快的优点。
RLS算法的关键是用二乘方的时间平均的最小化准则取代最小均方准则，并按时间迭代计算。具体来说，是要对初始时刻到当前时刻所有误差的平方进行平均并使其最小化，在按照这一准则确定FIR滤波器的权系数矢量，即所依据的准则是
(1)
其中式中，是期望响应，是L阶FIR滤波器的输出相应，即
(2)

2. 自适应滤波器系统辨识原理
如下图所示为自适应滤波器辨识未知系统的基本结构图，其中未知系统可以使FIR结构，也可以是IIR结构。在结构图中，x(n)为辨识系统而产生的输入信号，通常可以选择为白噪声，同时输入未知系统和自适应滤波器。d(n)为未知系统的输出信号，即自适应滤波器的参考信号。调整自适应滤波器的系数，使误差信号e(n)的均方误差达到最小，则自适应滤波器的输出y(n)近似等于未知系统的输出d(n)。可以证明，加性噪声v(n)的存在并不影响自适应滤波器最终收敛到最优维纳解。可以认为，具有相同输入和相似输出的两个系统，应该具有相似的特性。因此，可以采用自适应滤波器的特性或其单位脉冲响应来近似替代未知系统的特性或单位脉冲响应。

图1 自适应系统辨识原理图
3. Matlab仿真结果
下面给出使用RLS方法辨识未知系统的的Matlab仿真结果，并且修改其中的一些参数来比较不同情况下辨识系统的效果。

改变滤波器长度：

图2

图3

图4

图5

图6

图7

4. 相关matlab代码：

%采用RLS(递推最小二乘)方法辨识系统;
close all;
clc;
clear;

%辨识IIR类型系统;
NL=1000;
x=randn(1,NL); %输入信号
[N,Wc]=buttord([0.3,0.8],[0.4,07],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=30;
W0=zeros(1,M);
P0=5*eye(M);
lam=1;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(1);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差');
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');

%改变滤波器长度对辨识结果的影响;
%迭代次数为1000;
%遗忘因子为1;
%滤波器长度为10;
NL=1000;
x=randn(1,NL); %输入信号
[N,Wc]=buttord([0.3,0.8],[0.4,0.7],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=10;
W0=zeros(1,M);
P0=5*eye(M);
lam=1;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(2);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差');
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');

%迭代次数为1000;
%遗忘因子为1;
%滤波器长度为20;
NL=1000;
x=randn(1,NL); %输入信号
[N,Wc]=buttord([0.3,0.8],[0.4,0.7],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=20;
W0=zeros(1,M);
P0=5*eyeM);
lam=1;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(3);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差');
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');

%%改变遗忘因子对辨识结果的影响;
%迭代次数为1000;
%遗忘因子为0.97;
%滤波器长度为30;
NL=1000;
x=randn(1,NL); %输入信号
[N,Wc]=buttord([0.3,0.8],[0.4,0.7],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=30;
W0=zeros(1,M);
P0=5*eye(M);
lam=0.97;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(4);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差')
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');

%迭代次数为1000;
%遗忘因子为0.95;
%滤波器长度为30;
NL=1000;
x=randn(1,NL); %输入信号
[N,Wc]=buttord[0.3,0.8],[0.4,0.7],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=30;
W0=zeros(1,M);
P0=5*eye(M);
lam=0.95;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(5);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差');
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');

%更改迭代次数
%迭代次数为500;
%遗忘因子为1;
%滤波器长度为30;
NL=500;
x=randn(1,NL); %输入信号
[N,Wc]=buttord([0.3,0.8],[0.4,0.7],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=30;
W0=zeros(1,M);
P0=5*eye(M);
lam=1;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(6);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差');
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');

%更改迭代次数
%迭代次数为200;
%遗忘因子为1;
%滤波器长度为30;
NL=200;
x=randn(1,NL); %输入信号
[N,Wc]=buttord([0.3,0.8],[0.4,0.7],2,40);
[b,a]=butter(N,Wc,'stop'); % 需要辨识的IIR系统，此系统为一个巴特沃斯带通滤波器;
d=filter(b,a,x); % 参考信号
%初始化自适应滤波器系数;
M=30;
W0=zeros(1,M);
P0=5*eye(M);
lam=1;
S=initrls(W0,P0,lam);
[y,e,S]=adaptrls(x,d,S);
figure(7);
subplot(2,1,1);
plot(1:NL,[d;y;e]);
title('IIR系统的辨识');
xlabel('时间值');
ylabel('信号值');
legend('参考值','输出值','误差');
subplot(2,1,2);
[H1,W1]=freqz(b,a);
Hf1=abs(H1);
plot(W1,Hf1);
hold on;
[H2,W2]=freqz(S.coeffs,1);
Hf2=abs(H2);
plot(W2,Hf2,'r');
title('IIR系统和自适应滤波器的幅频响应曲线');
legend('IIR系统幅频响应曲线','自适应滤波器的幅频响应曲线');