Qurak

Make Code Talk

0%

自适应滤波器算法

Posted on Edited on In

目录

1. 自适应滤波器的原理

1.1 自适应滤波器介绍

在许多信号处理系统中，并没有信号的先验统计特性，在遇到信道均衡、回声消除以及其他因素之间的系统模型时，不能使用固定的参数处理。我们通过加入算法调节系统滤波系数，让系统具有自动调节的功能，我们称之为自适应滤波器。

自适应滤波器 主要特点：

  • 没有关于待提取信息的先验统计知识
  • 直接利用 利用观测数据在观测过程中不断递归更新
  • 是一个最优化的结果

自适应滤波器的模型结构如下图所示，该系统由两个输入端：输入信号 x(n)，期望信号 d(n)，两个信号进行比较等到一个误差信号 e(n)，通过自适应滤波器对滤波器参数进行调整，最终使得 e(n) 均方值最小。自适应滤波器可以利用前一时刻已得的滤波器参数产生的结果反馈到系统，调整参数，使得最终达到最优状态。

img

1.2 自适应滤波器模型分析

自适应滤波器模型如下图所示：

​ img

输入端矢量信号：x(n)

输出端矢量信号：y(n)

误差信号：e(n)

X(n)=[x(n),x(n−1),x(n−2)...,x(n−L)]TY(n)=∑k=0Lwk(n)x(n−k)=xT(n)w(n)e(n)=d(n)−y(n)=d(n)−wT(n)x(n) ξ(n)=E[e2(n)]=E[d2(n)]+wT(n)E[x(n)xT(n)]w(n)−2E[d(n)xT(n)]w(n)X(n)=[x(n),x(n-1),x(n-2)...,x(n-L)]^T\\ Y(n)=\sum^L_{k=0}w_k(n)x(n-k)=x^T(n)w(n)\\ e(n)=d(n)-y(n)=d(n)-w^T(n)x(n)\\ \ \\ \xi(n)=E[e^2(n)]=E[d^2(n)]+w^T(n)E[x(n)x^T(n)]w(n)-2E[d(n)x^T(n)]w(n)\\ \\

从这个式子可以看出，在输入信号和参考响应都是平稳随机信号的前提下，均方误差是权矢量的各个分量的 二次函数 。该函数图形是 L+2 维空间中的一个中间下凹的超抛物面，有唯一的最低点，称该曲面为 均方误差性能曲面。由此我们可以得到均方误差性能曲面的梯度：

∇=∂ξ∂w=2Rw−2P\nabla = \frac{\partial \xi}{\partial w}=2Rw-2P

梯度 ∇\nabla 等于0时，可求得最小均方误差对应的最佳权矢量或维纳解 W∗W^* ，解得 w∗=R−1Pw^*=R^{-1}P 。

为了得到这个解，我们还需要：

  • 事先知道 R、P 的值
  • 矩阵的逆计算量太大 O(n3)O(n^3)
  • 如果信号不是稳定的信号，则 R、P每次都不一样，需要重复计算

1.3 利用梯度下降算法和 LSM 算法

一般情况下，我们都采用梯度下降算法进行迭代搜索出最小值，得到的结果都是逼近维纳解而并非等于维纳解。利用最陡下降算法，沿着性能曲面负梯度最大的方向调整滤波器系数，以此来达到最低点。迭代公式：

w(n+1)=w(b)+μ(−∇ξ)w(n+1)=w(b)+\mu(-\nabla \xi)

μ\mu 为步长因子，μ\mu 的取值需要满足：（其中 λmax\lambda_{max} 表示输入信号自相关的最大特征值）

∣1−2μλmax∣<1⇒0<μ<1λmax|1-2\mu \lambda_{max}|<1\Rightarrow 0<\mu<\frac{1}{\lambda_{max}}

学过 线性代数 的同学知道，特征值可以用矩阵的迹来代替计算：

tr(R)=∑k=0Lλk>λmax⇒0<μ<1tr(R)<1λmaxtr(R)=\sum^L_{k=0}\lambda_{k}>\lambda_{max}\Rightarrow 0<\mu<\frac{1}{tr(R)}<\frac{1}{\lambda_{max}}

在梯度下降算法中，我们需要知道输入信号 x(n) 和期望信号 d(n) 的相关信息，当期望信号未知时，就无法确定他们的相关特性，必须对梯度向量进行估计。这里可以采用 LSM 自适应算法直接利用瞬态均方误差对瞬时抽头向量（滤波器系数）求梯度：

\begin{equation}\label{eq1} \hat\nabla \xi = \frac{\partial(e^2(n))}{\partial w(n)}=-2e(n)x(n) \end{equation}

所以传统 LMS 自适应滤波算法流程：

y(n)=wT(n)x(n)e(n)=d(n)−y(n)w(n+1)=w(n)+2μe(n)x(n)y(n)=w^T(n)x(n)\\ e(n)=d(n)-y(n)\\ w(n+1)=w(n)+2\mu e(n)x(n)\\

1.4 LMS 优缺点

LSM算法：

  • 优点：算法简单，易于实现，算法复杂度低，能够抑制旁瓣效应

  • 缺点：采用逐点更新，收敛速度慢；跟踪性能差，系统稳定性下降；LSM 要求不同时刻的输入向量 x 线性无关——LMS 的独立性假说，如果输入信号存在相关性，则会导致前一次迭代产生的梯度噪声传播到下一次迭代。

2 LMS 算法实现

LMS 自适应滤波器算法流程：

y(n)=wT(n)x(n)e(n)=d(n)−y(n)w(n+1)=w(n)+2μe(n)x(n)\\y(n)=w^T(n)x(n)\\ e(n)=d(n)-y(n)\\ w(n+1)=w(n)+2\mu e(n)x(n)\\

其中 w 为滤波系数，需要动态调节的；e 为代价函数（损失函数），用于梯度下降；u 为步长，需要自己设置。

2.1 Matlab 仿真

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
%------ LMS 滤波函数 ------
% 参数介绍：
%	xn:	输入信号
%	dn: 期望响应
%	M:	滤波器阶数
%	mu:	收敛因子（步长）
% 	itr:迭代次数
%	W:	滤波器系数矩阵
%	en:	误差序列
%	yn: 滤波器输出
function [yn, W, en] = LMSFilter(xn, dn, M, mu)
	itr = length(xn);
	en  = zeros(itr, 1);
	W   = zeros(M, itr);	% W 为一个 M x itr 的矩阵，初始化为每一次迭代初始值均为 0
	% 迭代计算
	for k = M : itr
		x = xn(k:-1:k-M+1)		% 滤波器前 M 个抽头输入
		y = W(: , k-1)‘.*x;		% 滤波器的输出
		en(k) = dn(k) - y;		% 第 k 次迭代的误差
		% 计算迭代式
		W(:,k) = W(:,k-1) + 2*mu * en(k) * x;
	end

	yn = inf*ones(size(xn));
	for k = M:length(xn)
		x = xn(k:-1:-M+1);
		yn(k) = W(:,end).'* x;  % 最终输出结果
	end
end

例程展示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
fs = 1000;	% 采样频率
f0 = 10;	% 信号频率
n = 500;	% 采样点数
t = (0:n-1)'/fs;
Signal0 = cos(2*pi*f0*t);
Signal1 = awgn(Signal0, 20);

figure;
subplot(211);plot(t,Signal0);title('原始信号');
subplot(212);plot(t,Signal1);title('加噪信号');

M = 30;
xn = Signal1;
dn = Signal0;

mu = 0.001;
[yn, W, en] = LMSFilter(xn, dn, M, mu);

figure;
ax1 = subplot(211);
plot(t, xn); grid on; ylabel('幅值')； xlabel('时间');
ylim([-1.5, 1.5]); title("LMS 滤波器输入信号");

ax2 = subplot(212);
plot(t, yn); grid on; ylabel('幅值'); xlabel('时间');
ylin([-1.5, 1.5]); title("LMS 滤波器输出信号");

figure; plot(en); grid on; title('误差');

2.2 C 语言代码

头文件 LMS.h

1
2
3
4
5
6
7
8
9
10
11
12
13
// ----- LMS.h ------
#ifndef _LMS_H_
#define _LMS_H_

#include <stdio.h>
#include <float.h>
#inlcude <math.h>
#include <string.h>

#define F_COUNT 1024
#define M 20

#endif

源文件 LMS.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
// ----- LMS.c ------
// xn	输入信号
// itr	迭代次数
// en	误差序列
// dn	期望响应
// M	滤波器阶数
// mu	收敛因子
// W	滤波器系数矩阵
// yn	实际输出
#include "LMS.h"
float *LMS_Filter(int itr, const float *xn, const float *dn, double mu, int length)
{
    static int i = 0;
    static int k = 0;
    static float y = 0.0;
    static float en[F_COUNT] = 0;		//	创建误差序列，并初始化 0
    static float W[M][F_COUNT] = 0;		//	创建系数矩阵，初始化为 0
    static float x[M] = 0;				
    static float yn[F_COUNT] = 0;
    
    for(k=M; k<=itr; k++)
    {
        // 滤波器 M 个抽头的输入：从 xn 第 k-1 个倒叙抽取 M 个样点的值放入 x
        // y 为滤波器输出： w 的第 k-2 列与 x 的积的和
        for(i=0; i<M; i++)
        {
            x[i]=xn[k-i-1];	// 	倒序抽取
            y += W[i][k-2] * x[i];
        }
        en[k-1] = dn[k-1] - y;	// 第 k 次迭代
        for(i=0; i<M; i++)
        {
            W[i][k-1] = W[i][k-2] + 2*mu*en[k-1]*x[i];
        }
        y = 0.0;
    }
    for(i=0; i<itr; i++)
    {        
        yn[i] = 0.0;
    }
    // 	求最优化时滤波器的输出序列
    for(k=M; k<= length; k++)
    {
        for(i=0; i<M; i++)
        {
            x[i] = xn[k-i-1];
            y += W[i][itr-1]*x[i];
        }
        yn[k-1] = y;
        y = 0.0;
    }
    
    return yn;
}

3 分块 LMS 算法（Block LMS）

3.1 BLMS 原理

相比于 LMS 算法中，每来一个采样点九调整一次滤波器权值，分块 LMS 算法则是每经过 k 个采样点才调整一次权值。其结构图如下图所示：

image-20210707045051184

输入数据序列经过转换器分成了若干个长度为 L 的数据块。然后将这些数据块一次一块的送入长度为 M 的滤波器中，在收集到每一块数据值后，进行滤波器抽头权值的更新，使得滤波器的自适应过程一块一块地进行，而不是像时域传统自适应滤波算法那样一个值一个值地进行

算法流程（积累 L 个采样点后更新滤波器系数）：

w(n+1)=w(n)+μX(n)e(n)w(n+2)=w(n+1)+μX(n+1)e(n+1)...w(n+L)=w(n+L−1)+μX(n+(L−1))e(n+(L−1))w(n+1) = w(n)+\mu X(n)e(n)\\ w(n+2) = w(n+1)+\mu X(n+1)e(n+1)\\ ...\\ w(n+L)=w(n+L-1)+\mu X(n+(L-1))e(n+(L-1))\\

所有式子相加得到：

w(n+L)=w(n)+μ∑i=0L−1X(n+i)e(n+i)w(n+L)=w(n)+\mu\sum_{i=0}^{L-1}X(n+i)e(n+i)

其中，μ\mu 为步长因子，更新一次 BLMS 相当于 LMS 更新了 L 次。例如采样 500 个点，L 选取 100，那么更新5 次 BLMS 相当于更新 500 次LMS，因此 BLMS 可以见小计算次数。

优点：相比于传统的 LMS 算法，运算量小一点

缺点：收敛速度与 LMS 算法相同

3.2 BLMS 算法实现

3.2.1 基于 Matlab 实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
% Block-Least-Mean-Square Algorithm
% d - 期望响应
% x - 输入信号
% mu - 步长因子
% M - 滤波器阶数
% L - 块长度

function [en, yn, w] = Block_LMS(d, x, mu, M, L)
	len_d = length(d);		% 输入信号
	k = floor(len_d / L);	% 分段，不足一段的部分舍弃
	y = zeros(len_d , 1);	% 初始化输出序列 y
	w = zeros(M, k+1);		% 滤波器系数矩阵: M 阶，递归 k 次
	e = zeros(len_d, 1);	% 初始化误差序列 e
	x_ = [zeros(M-1,1); x];
	
	for i = 1:k
		block_sum = 0;
		% 求每一个块的和 sum
		for j = (1+L*(k-1)):L*K
			X = x_(j:j+M-1);
			y(i) = w(:,k)'*x; % 滤波器输出
			
			e(i) = d(i) - y(i);
			block_sum = block_sum + X * e(i);
		end
		w(:k+1) = w(:,k) + mu * block_sum;
	end
	% 移动一步，使第 n 行对应第时间 n
	w = w(:,2:k+1);
end
3.2.2 基于 C 语言实现

头文件 BLMS.h

1
2
// ------ BLMS.h ------
// ...

源文件 BLMS.c

1
2
// ------ BLMS.c -------
// ...

4 BLMS 与 LMS 实际测试

在实际测试中其实可以发现 LMS 算法得到的结果：下面两个图分别表示 BLMS 与 LMS 算法在相同输入情况下的结果：

BLMS：

BLMS Algorithm

LMS：

LMS Algorithm

通过对比我们可以看出，BLMS 相对 LMS 而言收敛得更好，计算量更小。