rls_filter.m

来自「各种滤波器」· M 代码 · 共 66 行

function [y,w,NP]=RLS_Filter(x,d,N,delta,lambda,fashan)
% RLS_Filter 用于RLS滤波的函数,输入数据为x,期望数据为d
% 同时给出滤波器的阶数N,修正因子delta与遗忘因子lambda.
% 由于程序有两个算法,一个为数值发散的,一个是数值不发散的,所以
% fashan用来决定选用哪个算法.
% 默认值:N=10,delta=2,lambda=0.985,fashan=0;
% 函数返回滤波后的输出y和最后的滤波器系数w.
% 为了查看过程中矩阵P的范数,设置一个返回值NP,它与输入X同型,

% Argument check
if nargin<6
    fashan=0;
end
if nargin<5
    lambda=0.985;
end
if nargin<4
    delta=2;
end
if nargin<3
    N=10;
end
if nargin<2
    error('Input parameters must at least be 2 !');
end

% Parameters initialize
y=zeros(size(x));
w=zeros(N,1);
delaybuffer=zeros(N,1);
P=1/delta*eye(N);
NP=y;

if fashan==1
    % 下面的算法段来自《自适应滤波器原理》的P349
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for n=1:length(x)
        delaybuffer=[x(n);delaybuffer(1:end-1)];
        PAI=P*delaybuffer;
        k=PAI/(lambda+delaybuffer'*PAI);
        e=d(n)-w'*delaybuffer;
        w=w+k*e;
        g=k*delaybuffer';
        P=P/lambda-g/lambda*P;
        g=max(g(:));
        y(n)=w'*delaybuffer;
        NP(n)=norm(P);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    % 下面的算法段是对上面的算法段的改进版,改进了上面的P每次都要除以一个
    % 小于1的lambda带来的数值发散.下面的P(n)等价于上面的lambda^n*P(n)
    lambdan=lambda;
    for n=1:length(x)
        delaybuffer=[x(n);delaybuffer(1:end-1)];
        PAI=P*delaybuffer;
        k=PAI/(lambdan+delaybuffer'*PAI);
        e=d(n)-w'*delaybuffer;
        w=w+k*e;
        P=P-k*delaybuffer'*P;
        lambdan=lambdan*lambda;
        y(n)=w'*delaybuffer;
        NP(n)=norm(P);
    end
end