自适应滤波器（二）NLMS自适应滤波器

前一篇文章我们讲了LMS自适应滤波器，我们先回顾一下LMS算法流程：
yy(n)=wT(n)x(n)e(n)=d(n)−y(n)w(n+1)=w(n)+2μe(n)x(n)yy(n)=\boldsymbol{w}^{T}(n) \boldsymbol{x}(n) \\ \boldsymbol{e}(n)=d(n)-y(n) \\ \boldsymbol{w}(n+1)=\boldsymbol{w}(n)+2 \mu e(n) \boldsymbol{x}(n) yy(n)=wT(n)x(n)e(n)=d(n)−y(n)w(n+1)=w(n)+2μe(n)x(n)

影响LMS性能的因素，也就是最后一个公式的三个因素：

步长uuu，它是由我们事先指定
输入向量x(n)\mathbf{x}(n)x(n)
估计误差e(n)e(n)e(n)

如果x(n)\mathbf{x}(n)x(n)过大，那么2ue(n)x(n)2ue(n)\mathbf{x}(n)2ue(n)x(n)的结果中，x(n)\mathbf{x}(n)x(n)占的权重就太大了，而x(n)\mathbf{x}(n)x(n)是带噪信号，这样梯度噪声就被放大了。为了克服这个问题，可使用归一化LMS滤波器。在迭代时，对输入向量x(n)\mathbf{x}(n)x(n)*欧式范数（就是模值）*的平方进行归一化（Normalized LMS)。

归一化LMS滤波器是最小化干扰原理的一种表现形式，这个原理可以表述如下：

从一次迭代到下一次中，自适应滤波器的权向量应当以最小方式改变，而且受到更新的滤波器输出所施加的约束。

用w^(n)\hat{\mathbf{w}}(n)w^(n)表示第n次迭代滤波器的权向量，w^(n+1)\hat{\mathbf{w}}(n+1)w^(n+1)表示第n+1次迭代滤波器的权向量，那么NLMS设计准则可表述为约束优化问题：给定输入向量x(n)\mathbf{x}(n)x(n)和目标响应d(n)d(n)d(n)，确定更新抽头向量w^(n+1)\hat{\mathbf{w}}(n+1)w^(n+1)，以使如下增量
δw^(n+1)=w^(n+1)−w^(n)\delta \hat{\mathbf{w}}(n+1)=\hat{\mathbf{w}}(n+1)-\hat{\mathbf{w}}(n) δw^(n+1)=w^(n+1)−w^(n)
的欧式范数最小化，并受制于以下约束条件
w^H(n+1)x(n)=d(n)\hat{\mathbf{w}}^{\mathrm{H}}(n+1) \mathbf{x}(n)=d(n) w^H(n+1)x(n)=d(n)
使用拉格朗日乘子法来解决这个约束问题，那么代价函数为：
J(n)=∥δw^(n+1)∥2+Re⁡[λ∗(d(n)−w^H(n+1)x(n))]J(n)=\|\delta \hat{\mathbf{w}}(n+1)\|^{2}+\operatorname{Re}\left[\lambda^{*}\left(d(n)-\hat{\mathbf{w}}^{\mathrm{H}}(n+1) \mathbf{x}(n)\right)\right] J(n)=∥δw^(n+1)∥2+Re[λ∗(d(n)−w^H(n+1)x(n))]
其中，λ\lambdaλ为复数拉格朗日乘子，∗*∗表示复共轭，Re[⋅]Re[\cdot]Re[⋅]表示取实部运算，约束对代价函数的贡献是实值的；∥δw^(n+1)∥2\|\delta \hat{\mathbf{w}}(n+1)\|^{2}∥δw^(n+1)∥2表示欧式范数的平方运算，其结果也是实数。因此，代价函数J(n)J(n)J(n)是实值的二次函数，且可表示为：
J(n)=(w^(n+1)−w^(n))H(w^(n+1)−w^(n))+Re⁡[λ∗(d(n)−w^H(n+1)x(n))]J(n)=(\hat{\mathbf{w}}(n+1)-\hat{\mathbf{w}}(n))^{\mathrm{H}}(\hat{\mathbf{w}}(n+1)-\hat{\mathbf{w}}(n))+\operatorname{Re}\left[\lambda^{*}\left(d(n)-\hat{\mathbf{w}}^{\mathrm{H}}(n+1) \mathbf{x}(n)\right)\right] J(n)=(w^(n+1)−w^(n))H(w^(n+1)−w^(n))+Re[λ∗(d(n)−w^H(n+1)x(n))]
采用如下步骤寻找最小的最优更新权向量：

代价函数J(n)J(n)J(n)对w^(n+1)\hat{\mathbf{w}}(n+1)w^(n+1)求导，可得：

∂J(n)∂w^H(n+1)=2(w^(n+1)−w^(n))−λ∗x(n)\frac{\partial J(n)}{\partial \hat{\mathbf{w}}^{\mathrm{H}}(n+1)}=2(\hat{\mathbf{w}}(n+1)-\hat{\mathbf{w}}(n))-\lambda^{*} \mathbf{x}(n) ∂w^H(n+1)∂J(n)=2(w^(n+1)−w^(n))−λ∗x(n)

令其为零，即得最优解为：
w^(n+1)=w^(n)+12λ∗x(n)\hat{\mathbf{w}}(n+1)=\hat{\mathbf{w}}(n)+\frac{1}{2} \lambda^{*} \mathbf{x}(n) w^(n+1)=w^(n)+21λ∗x(n)

将上式带入约束条件，求解未知乘子λ\lambdaλ：

d(n)=w^H(n+1)x(n)=(w^(n)+12λ∗x(n))Hx(n)=w^H(n)x(n)+12λuH(n)x(n)=w^H(n)x(n)+12λ∥x(n)∥2\begin{aligned} d(n) &=\hat{\mathbf{w}}^{\mathrm{H}}(n+1) \mathbf{x}(n) \\ &=\left(\hat{\mathbf{w}}(n)+\frac{1}{2} \lambda^{*} \mathbf{x}(n)\right)^{\mathrm{H}} \mathbf{x}(n) \\ &=\hat{\mathbf{w}}^{\mathrm{H}}(n) \mathbf{x}(n)+\frac{1}{2} \lambda \mathbf{u}^{\mathrm{H}}(n) \mathbf{x}(n) \\ &=\hat{\mathbf{w}}^{\mathrm{H}}(n) \mathbf{x}(n)+\frac{1}{2} \lambda\|\mathbf{x}(n)\|^{2} \end{aligned} d(n)=w^H(n+1)x(n)=(w^(n)+21λ∗x(n))Hx(n)=w^H(n)x(n)+21λuH(n)x(n)=w^H(n)x(n)+21λ∥x(n)∥2

可求得：
λ=2e(n)∥u(n)∥2\lambda=\frac{2 e(n)}{\|\mathbf{u}(n)\|^{2}} λ=∥u(n)∥22e(n)
其中，
e(n)=d(n)−w^H(n)x(n)e(n)=d(n)-\hat{\mathbf{w}}^{\mathrm{H}}(n) \mathbf{x}(n) e(n)=d(n)−w^H(n)x(n)
是误差信号。

结合前两步的结果，可得：

δw^(n+1)=w^(n+1)−w^(n)=1∥x(n)∥2x(n)e∗(n)\begin{aligned} \delta \hat{\mathbf{w}}(n+1) &=\hat{\mathbf{w}}(n+1)-\hat{\mathbf{w}}(n) \\ &=\frac{1}{\|\mathbf{x}(n)\|^{2}} \mathbf{x}(n) e^{*}(n) \end{aligned} δw^(n+1)=w^(n+1)−w^(n)=∥x(n)∥21x(n)e∗(n)

为了对一次迭代到下一次迭代抽头权向量的增量变化进行控制而不改变向量的方向，引入一个正的实数标度因子uuu，该增量可以写为：
δw^(n+1)=w^(n+1)−w^(n)=μ~∥x(n)∥2x(n)e∗(n)\begin{aligned} \delta \hat{\mathbf{w}}(n+1) &=\hat{\mathbf{w}}(n+1)-\hat{\mathbf{w}}(n) \\ &=\frac{\widetilde{\mu}}{\|\mathbf{x}(n)\|^{2}} \mathbf{x}(n) e^{*}(n) \end{aligned} δw^(n+1)=w^(n+1)−w^(n)=∥x(n)∥2μx(n)e∗(n)
等价的，我们可以写出：
w^(n+1)=w^(n)+μ∥x(n)∥2x(n)e∗(n)\hat{\mathbf{w}}(n+1)=\hat{\mathbf{w}}(n)+\frac{\mu}{\|\mathbf{x}(n)\|^{2}} \mathbf{x}(n) e^{*}(n) w^(n+1)=w^(n)+∥x(n)∥2μx(n)e∗(n)
这个公式就是归一化LMS算法抽头权向量的递归公式，为什么叫归一化呢？因为公式中对输入向量$ \mathbf{x}(n)$欧式范数的平方就行了归一化。

当输入向量$\mathbf{x}(n) 较小时，较小时，较小时，\frac{\mu}{|\mathbf{x}(n)|^{2}}$的值过小，可能导致数值计算困难的情况，为了克服这个情况，将上面的表达式改为：
w^(n+1)=w^(n)+μ~δ+∥x(n)∥2x(n)e∗(n)\hat{\mathbf{w}}(n+1)=\hat{\mathbf{w}}(n)+\frac{\widetilde{\mu}}{\delta+\|\mathbf{x}(n)\|^{2}} \mathbf{x}(n) e^{*}(n) w^(n+1)=w^(n)+δ+∥x(n)∥2μx(n)e∗(n)
其中，δ>0\delta>0δ>0。

我们总结NLMS算法的步骤如下：

如果有抽头权向量w^(n)\hat{\mathbf{w}}(n)w^(n)的先验知识，则用它来做初值，否则令w^(n)\hat{\mathbf{w}}(n)w^(n)为0

计算输出：y(n)=wT(n)x(n)y(n)=\boldsymbol{w}^{T}(n) \boldsymbol{x}(n)y(n)=wT(n)x(n)

计算误差：e(n)=d(n)−y(n){e}(n)=d(n)-y(n)e(n)=d(n)−y(n)

权重更新：w^(n+1)=w^(n)+μ~δ+∥x(n)∥2x(n)e∗(n)\hat{\mathbf{w}}(n+1)=\hat{\mathbf{w}}(n)+\frac{\widetilde{\mu}}{\delta+\|\mathbf{x}(n)\|^{2}} \mathbf{x}(n) e^{*}(n)w^(n+1)=w^(n)+δ+∥x(n)∥2μx(n)e∗(n)

% 输入参数:
%   xn   输入的信号，列向量
%   dn   所期望的响应
%   M    滤波器的阶数
%   mu   收敛因子(步长)
% 输出参数:
%   W    滤波器系数矩阵
%   en   误差序列
%   yn   滤波器输出
function [yn, W, en]=nlmsFunc(xn, dn, M, mu, delta)
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) + mu*en(k)*x/(delta+x.' * x);
endyn = inf * ones(size(xn)); % 初值为无穷大是绘图使用，无穷大处不会绘图
for k = M:length(xn)x = xn(k:-1:k-M+1);yn(k) = W(:,end).'* x;  % 最终输出结果
end

%% 产生测试信号
fs = 1;
f0 = 0.02;
n = 1000;
t = (0:n-1)'/fs;
xs = cos(2*pi*f0*t);
ws = awgn(xs, 20, 'measured');M  = 20 ;   % 滤波器的阶数
xn = ws;
dn = xs;
% rho_max = max(eig(ws*ws.'));   % 输入信号相关矩阵的最大特征值
% mu = (1/rho_max) ;    % 收敛因子 0 < mu < 1/rho
mu1 = 0.01;
mu2 = 1;
delta = 1e-3;
[yn,W,en] = lmsFunc(xn,dn,M,mu1);
[yn2,W2,en2] = nlmsFunc(xn, dn, M, mu2, delta);

欢迎关注微信公众号：Quant_Times

欢迎大家学习我的课程：
System Generator & HLS数字信号处理教程

自适应滤波器（二）NLMS自适应滤波器相关推荐

自适应滤波：最小均方误差滤波器（LMS、NLMS）
作者:桂. 时间:2017-04-02 08:08:31 链接:http://www.cnblogs.com/xingshansi/p/6658203.html 声明:欢迎被转载,不过记得注明出处哦 ...
[转载+原创]Emgu CV on C# （五） —— Emgu CV on 局部自适应阈值二值化
局部自适应阈值二值化相对全局阈值二值化,自然就有局部自适应阈值二值化,本文利用Emgu CV实现局部自适应阈值二值化算法,并通过调节block大小,实现图像的边缘检测. 一.理论概述(转载自< ...
二维高斯滤波器（gauss filter）的实现
我们以一个二维矩阵表示二元高斯滤波器,显然此二维矩阵的具体形式仅于其形状(shape)有关: def gauss_filter(kernel_shape): 为实现二维高斯滤波器,需要首先定义二元高斯 ...
LabView学习笔记（二）：滤波器实验
Labview学习笔记: LabView学习笔记(一):基础介绍 LabView学习笔记(二):滤波器实验 LabView学习笔记(三):基本控件 LabView学习笔记(四):动态数据类型 LabV ...
python对参数二值化处理_OpenCV自适应阀值二值化表格检测方法(python版)
OCR主要分为三个步骤:检测.分割.文字识别.其中文字识别无论是英文还是中文相对比较成熟.只要检测到位,标准的印刷体识别率还是非常高的. 文书OCR检测主要有文字检测和表格检测.文本段落基于行的检测通 ...
空时二维自适应处理-相控阵雷达二维杂波谱分布仿真Matlab
空时二维自适应处理-相控阵雷达二维杂波谱分布仿真非均匀性杂波回波相控阵雷达二维杂波谱分布相控阵雷达二维杂波谱分布仿真非均匀性杂波回波在机载雷达杂波仿真时,需要综合考虑杂波单元的反射率模型 ...
组合导航系列文章（十二）：滤波器基本原理
我觉得不错的地方 3. 滤波器估计原理组合导航中,先验是imu解算的值,观测是gps等传感器给出的值,融合的目的是找到概率最大的那个值.上面介绍的三种方法都是对先验和观测的融合,由于极大似然要求先验 ...
免费分享thinkphp框架开发周易八字起名网宝宝起名在线下单网站源码自适应可二开
宝宝起名/八字起名/周易取名/周易八字起名平台网站/在线付费起名源码,thinkphp框架开发周易八字起名网宝宝起名在线下单网站源码自适应可二开,PHP权威起名策划机构平台源码,Thinkphp3.2 ...
MAX30102脉搏血氧仪和心率传感器（二）FIR滤波器
文章目录前言一.修正上一章产生的错误二.FIR滤波器设计 1.对采集的信号进行频谱分析 2.滤波器设计 3.滤波器仿真三.ARM_MATH库实现(以STM32为例) 实际效果测试滤波前滤波 ...
中点和中值滤波的区别_组合导航系列文章（十二）：滤波器基本原理
<组合导航系列文章>是<从零开始做自动驾驶定位>系列的第二阶段,从本阶段开始,文章在<泡泡机器人>公众号上首发,知乎用来备份和以后可能出现的必要更正. 泡泡机器人文 ...

自适应滤波器（二）NLMS自适应滤波器

自适应滤波器（二）NLMS自适应滤波器相关推荐

最新文章

热门文章