浅谈自适应滤波器---（快速RLS算法）

在上一篇博客中（浅谈自适应滤波器）我给大家介绍了关于自适应滤波器的一些入门级的知识，并分析了常规RLS算法单次迭代的计算量级为O[N2]，当阶数N增大时相应的计算量显著增大，为了将计算量级降低到O[N]，人们提出了很多快速的RLS算法，其中快速横向递归最小二乘算法(FTRLS)最具有吸引力，得到了广泛的应用。
本次将为大家介绍该算法的由来并给出相应的算法，在推导FTRLS算法时，需要同时求解前向和后向线性预测问题，以及另外两个横向滤波器：联合过程估计器和辅助滤波器；实际上快速横向RLS算法可以视为这四个并行工作的滤波器，它们之间互相交换参量，相互作用最后得到RLS的快速实现方法。

该图显示了这四个滤波器的关系。

---

向前预测关系式

对于一个N阶的预测器，其瞬时后验前向预测误差可以表示为

后验和先验前向预测误差之间的关系有下式给出（所谓的先验是指利用k-1时刻的滤波系数矢量估计k时刻的输出，而后验是指利用k时刻的滤波系数矢量估计k时刻的输出）

其中gamma(k-1,N)称为k-1时刻N阶的转换因子（表示公式中的希腊字母“gamma”，下同）。
最小加权最小二乘误差的时间更新方程如下

还有转换因子的更新方程

接下来就是前向预测抽头系数向量更新方程

其中Phi(k,N+1)的更新方程为

后向预测关系式

同样后验和先验的后向预测误差关系表达式为

最小加权最小二乘误差的时间更新方程如下

后向预测抽头系数向量更新方程

其中Phi(k,N)的更新方程为

则Phi(k,N+1)中的最后一个元素可以表示为

最后就是关于gamma(k,N)更新方程

联合过程估计

该过程实现了输入信号与期望信号的联合估计，其先验误差可写成

其先验误差与后验误差的关系为

联合过程估计其抽头系数的时间更新关系式为

至此算法所需的计算表达式已全部给出，下面给出的是经过整理完整的算法描述。

---

Matlab仿真

首先给出m代码

%FTF快速横向滤波算法
clc;
clear all;
close all;
%************************生成仿真信号**************************************
Fs = 10000;                                                     %设置采样频率
t = 0:1/Fs:3.5;  
t = t';
Size_t = size(t,1);
F1 = 2;
F2 = 10;
F3 = 20;
F4 = 1000;
Signal = sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t); %生成信号
noise_amp = 1;                                           %定义噪声的标准差
noise1 = noise_amp*randn(Size_t,1);                        %生成高斯白噪声
noise2 = noise_amp*randn(Size_t,1);
noise3 = 5*sin(2*pi*F4*t);

noise = noise2;
Signal_noise = Signal + 0.2*noise;                           %加入高斯白噪声
Signal_noise(2:end) = Signal_noise(2:end) + 0.15*noise(1:end-1);
Signal_noise(3:end) = Signal_noise(3:end) + 0.1*noise(1:end-2);

subplot(2,1,1);
plot(t,Signal);
title('原始信号');
subplot(2,1,2);
plot(t,Signal_noise);
title('加入干扰噪声的信号');
%*************************************************************************
N = 3;                     %定义FIR滤波器阶数
Signal_Len = Size_t;       %定义信号数据的个数
lambda = 1;                %定义遗忘因子
delta = 0.01;    
y_out = zeros(Signal_Len,1);
Eta_out = zeros(Signal_Len,1);
w_out = zeros(Signal_Len,N);
for i = 1:Signal_Len
    if i==1
        w_f_last = zeros(N,1);
        w_b_last = zeros(N,1);
        w_last = zeros(N,1);
        Phi_last = zeros(N,1);
        gamma_last = 1;
        xi_f_last = delta;
        xi_b_last = delta;
        x_N_1 = zeros(N+1,1);
    end
    %输入数据
    if i<= N
        x_N_1(1:i) = noise(i:-1:1,1);
    else        
        x_N_1 = noise(i:-1:(i-N),1);
    end    
    d = Signal_noise(i);
    %算法主体
    e_f = x_N_1' * [1;-w_f_last];                                          %(1)
    epsilon_f = e_f * gamma_last;                                          %(2)
    xi_f = lambda * xi_f_last + e_f * epsilon_f;                           %(3)
    w_f = w_f_last + Phi_last * epsilon_f;                                 %(4)
    Phi_N_1 = [0;Phi_last] + e_f/(lambda * xi_f_last)*[1;-w_f_last];       %(5)
    Phi_N_1_N_1 = Phi_N_1(end);
    gamma_N_1 = (lambda * xi_f_last)/xi_f * gamma_last;                    %(6)
    e_b = lambda * xi_b_last * Phi_N_1_N_1;                                %(7)
    gamma = 1/(1/gamma_N_1 - Phi_N_1_N_1*e_b);                             %(8)
    epsilon_b = e_b * gamma;                                               %(9)
    xi_b = lambda * xi_b_last + epsilon_b * e_b;                           %(10)
    Phi = Phi_N_1 - Phi_N_1_N_1 * [-w_b_last;1];                           %(11)
    Phi = Phi(1:end-1);                                                    
    w_b = w_b_last + Phi * epsilon_b;                                      %(12)
    x = x_N_1(1:N);
    y = w_last'* x;                                                  
    e = d - y;                                                             %(13)
    epsilon = e * gamma;                                                   %(14)
    w = w_last + Phi * epsilon;                                            %(15)
    %变量更替
    xi_f_last = xi_f;
    w_f_last = w_f;
    gamma_last = gamma;
    xi_b_last = xi_b;
    Phi_last = Phi;
    w_b_last = w_b;
    w_last = w;
    %滤波结果存储
    y_out(i) = y;
    Eta_out(i) = e;
    w_out(i,:) = w';
end
figure;
subplot(2,1,1);
plot(y_out);
title('滤波器输出');
subplot(2,1,2);
plot(Eta_out);
title('输出误差');

figure;
plot(t(1:Signal_Len),w_out(:,1),'r',t(1:Signal_Len),w_out(:,fix(N/2)+1),'b',t(1:Signal_Len),w_out(:,N),'y');
title('自适应滤波器系数');

运行一下得到如下的效果图



对照我在上一篇博客里给出的常规RLS结果，可以发现一个直观的结果就是滤波的效果比常规的稍微差一些，但比LMS的还是要好很多，同时单次迭代的计算量和LMS的相当，事实上可以通过改变代码里面的一些参数滤波的效果又会得到进一步的改善，因为此快速算法是在假定运算精度无穷高的情况下得到的，而在Matlab仿真时计算的精度是有限的（即使是double型，也是有限的）因此会与理论的有偏差，这些只有通过实验才能确切的感受到。

---

到此为止，我们就可以利用上面我给大家介绍的算法做一些比较高大上的项目了，我最近一直在思考怎么来把这些理论运用到实际的工程当中，想来想去还真的想出了一个比较靠谱的点子，这要得益于本人本科时在科创方面的豪投入，我手里正好有一块闲置很久的TMS320F2812开发板，最关键的是上面有现成的ADC输入和DAC输出接口，我准备利用它来做一个能够消除声学回声的扩音器系统。这实际上就是典型应用中的信号增强部分，我在后面的博客中会对这个项目进行跟进，向大家展示自适应滤波器的强大魅力。