在很多信号处理系统中,并没有信号的先验统计特性,不能使用某一固定参数的滤波器来处理,比如信道均衡、回声消除以及其他因素之间的系统模型等,均采用了调整系数的滤波器,称为自适应滤波器。这样的滤波器结合了允许滤波器系数适应于信号统计特性的算法。
自适应滤波器的原理框图如下图所示,输入信号x(n) 通过参数可调数字滤波器后产生输出信号 y(n),将其与期望信号d(n)进行比较,形成误差信号e(n), 通过自适应算法对滤波器参数进行调整,最终使 e(n)的均方值最小。自适应滤波可以利用前一时刻已得的滤波器参数的结果,自动调节当前时刻的滤波器参数,以适应信号和噪声未知的或随时间变化的统计特性,从而实现最优滤波。自适应滤波器实质上就是一种能调节自身传输特性以达到最优的维纳滤波器。自适应滤波器不需要关于输入信号的先验知识,计算量小,特别适用于实时处理。维纳滤波器参数是固定的,适合于平稳随机信号。卡尔曼滤波器参数是时变的,适合于非平稳随机信号。然而,只有在信号和噪声的统计特性先验已知的情况下,这两种滤波技术才能获得最优滤波。在实际应用中,常常无法得到信号和噪声统计特性的先验知识。在这种情况下,自适应滤波技术能够获得极佳的滤波性能,因而具有很好的应用价值。
image-20210309162643639
自适应滤波器的详细模型如下图所示。
image-20210310142959048
输入信号矢量:
输出为:
自适应线性组合器的L+1个权系数构成一个权系数矢量,称为权矢量,用表示,即
因此可以表示为:
误差信号为:
自适应线性组合器按照误差信号均方值最小的准则,即:
输入信号的自相关矩阵为:
期望信号与输入信号的互相关矩阵为:
则均方误差的简单表示形式为:
从该式可看出,在输入信号和参考响应都是平稳随机信号的前提下,均方误差是权矢量的各分量的二次函数。该函数图形是L+2维空间中一个中间下凹的超抛物面,有唯一的最低点,该曲面称为均方误差性能曲面,简称性能曲面。
均方误差性能曲面的梯度:
令梯度 等于零,可求得最小均方误差对应的最佳权矢量或维纳解 ,解得 }=\boldsymbol{R}^{-1} \boldsymbol{P}w∗=R−1P。
虽然维纳解的表达式我们知道了,但仍然有几个问题:
2.3 梯度下降法
一般情况下,我们使用递归的方法来寻找多变量函数的最小值,其性能指标就是MSE(Mean Square Error),它是滤波器系数的二次函数,因此该函数具有唯一的最小值。一般是采用梯度下降的方法来进行迭代搜索出最小值,梯度下降又分为梯度下降、随机梯度下降和批量梯度下降。
使用迭代搜索的方式一般都只能逼近维纳解,并不等同于维纳解。
下面我们先来看梯度下降法,再来看下前面的公式:(梯度下降的原理可参考我的另一篇文章:基于梯度下降算法的线性回归拟合(附python/matlab/julia代码))
误差信号为:
均方误差为:
利用最陡下降算法,沿着性能曲面最速下降方向(负梯度方向)调整滤波器强权向量,搜索性能曲面的最小点,计算权向量的迭代公式为:
其中为步长因子,的取值需要满足下式,其中表示输入信号自相关矩阵的最大特征值。
由于计算特征值比较复杂,有时为了避免计算特征值,可采用计算矩阵迹的方法,因为自相关矩阵是正定的,因此有:
这里,是的迹,它可以用输入信号的取样值进行估计,即:
因此有:
在梯度下降算法中,为获得系统的最佳维纳解,需要知道输入信号和期望信号的相关信息,当期望信号未知时,就无法确定它们的相关特性,必须对梯度向量进行估计。
LMS自适应算法直接利用瞬态均方误差对瞬时抽头向量(滤波器系数)求梯度:
由此可得传统LMS自适应滤波算法流程如下:
下面我们来看一个实际的仿真:
close all;clear all;clc;
%% 产生测试信号
fs = 1;
f0 = 0.02;
n = 1000;
t = (0:n-1)'/fs;
xs = cos(2*pi*f0*t);
ws = awgn(xs, 30, 'measured');
% figure;
% subplot(211);plot(t, xs);title('原始信号');
% subplot(212);plot(t, ws);title('加噪信号');
M = 20 ; % 滤波器的阶数
xn = ws;
dn = xs;
% rho_max = max(eig(ws*ws.')); % 输入信号相关矩阵的最大特征值
% mu = (1/rho_max) ; % 收敛因子 0 < mu < 1/rho
mu = 0.001;
[yn,W,en] = lmsFunc(xn,dn,M,mu);
figure;
ax1 = subplot(211);
plot(t,ws);grid on;ylabel('幅值');xlabel('时间');
ylim([-1.5 1.5]);title('LMS滤波器输入信号');
ax2 = subplot(212);
plot(t,yn);grid on;ylabel('幅值');xlabel('时间');title('LMS滤波器输出信号');
ylim([-1.5 1.5]);linkaxes([ax1, ax2],'xy');
figure;plot(en);grid on;title('误差');
其中LMS滤波器函数如下:
% 输入参数:
% xn 输入的信号
% dn 所期望的响应
% M 滤波器的阶数
% mu 收敛因子(步长)
% 输出参数:
% W 滤波器系数矩阵
% en 误差序列
% yn 滤波器输出
function [yn, W, en]=lmsFunc(xn, dn, M, mu)
itr = length(xn);
en = zeros(itr,1);
W = zeros(M,itr); % 每一列代表-次迭代,初始为0
% 迭代计算
for k = M:itr % 第k次迭代
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:k-M+1);
yn(k) = W(:,end).'* x; % 最终输出结果
end
image-20210317123716776
这里可能有同学对期望信号会有疑惑,期望信号都是已知的了,还滤波干嘛?LMS滤波器还有什么用?
;
LMS算法的优缺点:
正是由于LMS算法的缺陷,后面才有了NLMS、RLS等算法,我们会在后面的文章中一一讲到。
附:上述仿真的Python代码如下:
# This is a sample Python script.
# Press Shift+F10 to execute it or replace it with your code.
# Press Double Shift to search everywhere for classes, files, tool windows, actions, and settings.
import numpy as np
import matplotlib.pyplot as plt
import scipy.io as scio
# 信号加噪
def awgn(x, snr):
snr = 10 ** (snr / 10.0)
xpower = np.sum(np.abs(x) ** 2) / len(x)
npower = xpower / snr
if type(x[0]) != np.complex128:
return x + np.random.randn(len(x)) * np.sqrt(npower)
else:
return x + np.random.randn(len(x)) * np.sqrt(npower / 2) + 1j * np.random.randn(len(x)) * np.sqrt(npower / 2)
def lmsFunc(xn, dn, M, mu):
itr = len(xn)
en = np.zeros((itr, 1))
W = np.zeros((M, itr))
for k in range(M, itr):
if k==20:
x = xn[k-1::-1]
else:
x = xn[k-1:k-M-1:-1]
try:
y = np.dot(W[:, k - 2], x)
print(y)
except:
pass
en[k-1] = dn[k-1] - y
W[:, k-1] = W[:, k - 2] + 2 * mu * en[k-1] * x
yn = np.ones(xn.shape) * np.nan
for k in range(M, len(xn) ):
if k == 20:
x = xn[k - 1::-1]
else:
x = xn[k - 1:k - M - 1:-1]
yn[k] = np.dot(W[:, -2], x)
return yn, W, en
if __name__ == '__main__':
fs = 1
f0 = 0.02
n = 1000
t = np.arange(n)/fs
xs = np.cos(2*np.pi*f0*t)
ws = awgn(xs, 20)
# data1 = scio.loadmat('xs.mat')
# data2 = scio.loadmat('ws.mat')
# xs = data1['xs'].flatten()
# ws = data2['ws'].flatten()
M = 20
xn = ws
dn = xs
mu = 0.001
yn, W, en = lmsFunc(xn, dn, M, mu)
plt.figure()
plt.subplot(211)
plt.plot(t, ws)
plt.subplot(212)
plt.plot(t, yn)
plt.figure()
plt.plot(en)
plt.show()