自适应滤波算法原理及其应用.doc

参考文献
自适应滤波算法原理与应用
经典的滤波算法包括，维纳滤波，卡尔曼滤波，自适应滤波。维纳滤波与卡尔曼滤波能够满足一些工程问题的需求，得到较好的滤波效果。但是他们也存在局限性，对于维纳滤波来说，需要得到足够多的数据样本时，才能获得较为准确的自相关函数估计值，一旦系统设计完毕，滤波器的长度就不能再改变，这难以满足信号处理的实时性要求；对于卡尔曼滤波，需要提前对信号的噪声功率进行估计，参数估计的准确性直接影响到滤波的效果。在实际的信号处理中，如果系统参数能够随着输入信号的变化进行自动调整，不需要提前估计信号与噪声的参数，实现对信号的自适应滤波，这样的系统就是自适应滤波系统。
1.基本自适应滤波算法
自适应滤波算法的基本思想是根据输入信号的特性自适应调整滤波器的系数，实现最优滤波。

图1 自适应滤波结构框图
若自适应滤波的阶数为，滤波器系数为，输入信号序列为，则输出为：
（ 1）
（ 2）
其中为期望信号，为误差信号。
（ 3）
（ 4）
则滤波器的输出可以写成矩阵形式：
（ 5）
（ 6）
定义代价函数：
（ 7）
当使上式中的代价函数取到最小值时，认为实现最优滤波，这样的自适应滤波成为最小均方自适应滤波（LMS）。
对于最小均方自适应滤波，需要确定使得均方误差最小的滤波器系数，一般使用梯度下降法求解这类问题。滤波器系数向量的迭代公式为：
（ 8）
式中，为步长因子，为代价函数的梯度。
（ 9）
因为瞬时梯度为真实梯度值的无偏估计，实际应用中可使用瞬时梯度代替真实梯度，即有：
（ 10）
（ 11）
通过逐步迭代，即可得到最优的滤波器系数，实现对输入信号的自适应滤波。
2.自适应滤波的工程应用
为了比较不同滤波算法的滤波效果，这里仍然采用前面用到的二维圆周运动轨迹追踪的问题作为工程背景。自适应滤波算法的程序设计思路如图 2所示。
图 2 自适应滤波算法流程图
迭代步长时，得到的滤波结果为：

图 3 X方向自适应滤波结果-基本自适应滤波

图 4 Y方向自适应滤波结果-基本自适应滤波
从X与Y方向上的位移变化曲线与方差变化曲线上可以看出，滤波结果出现了发现，最终得到的结果并没有达到最优解。分析其原因，可能是迭代步长太大，将迭代步长减小之后，取，得到较为理想的滤波结果，示于图 5和 6.

图 5 X方向自适应滤波结果-基本自适应滤波

图 6 Y方向自适应滤波结果-基本自适应滤波
可以看出，减小步长因子之后，两个方向上的滤波轨迹与期望的轨迹之间的误差明显减小，证明了自适应滤波的有效性。
3.自适应滤波的收敛性分析
在上一节的讨论中，迭代步长选择对于算法的收敛性具有决定性作用，步长值的微小改变即可对算法的收敛效果产生明显影响，因此如何确定合适的步长值是自适应滤波算法中重要的内容。
（ 12）
（ 13）
系统的最小均方误差最小时，有：
则下式成立： （ 14）
对于滤波器系数的迭代过程，有：
（ 15）
对自相关矩阵进行分解，即： （ 16）
则相邻两次迭代过程的滤波器系数之间满足关系式：
（ 17）
（ 18）
当迭代次数为无穷大时，理论上可以实现最优滤波，即迭代步长应该满足：
（ 19）
从而有：
（ 20）
式 20即为确保算法收敛迭代步长应满足的条件。
得到步长的收敛性条件，即可在满足要求的范围内调整步长因子，选择最佳的步长，在确保算法收敛的前提下，提高收敛速度。对于二维轨迹追踪问题，取步长因子为，得到的滤波结果如图 7至 9所示。

图 7 X方向自适应滤波结果-基本自适应滤波

图 8 Y方向自适应滤波结果-基本自适应滤波

图 9 二维圆周运动轨迹滤波结果-基本自适应滤波
从X方向，Y方向上的滤波结果可以看出，滤波轨迹在起初的一段时间内与期望轨迹存在较大的误差，但随着迭代次数增加，两者的误差逐渐减小，最终得到误差的最小值。二维轨迹图上也能得到类似的结论。
4.变步长自适应滤波
在满足收敛性条件的要求下选择迭代步长，可以确保最终得到收敛的结果，但是这一步长在整个过程中是固定的。然而，更为理想的情况是在滤波的初始阶段，误差值很大时，迭代步长可以取较大的值，以取得较快的收敛速度，随着误差减小，逐渐接近最优目标时，迭代步长也相应减小，从而得到较好的收敛精度，这就是变步长自适应滤波算法。
变步长的自适应滤波算法已经有了较长时间的发展，前人发展了很多有效的变步长算法，这里仅选择两种常用的方法。
（1）归一化变步长自适应滤波算法

其中，为常数，且满足。
归一化的变步长滤波算法使用输入信号的能量对步长因子进行归一化，确保其取到合适的值。
（2）Sigmod函数变步长自适应滤波算法

其中，为常数，且满足。
Sigmod函数使用滤波器的输出误差对迭代步长进行控制，从表达式中可以看出，误差较大时，步长因子的值较大，误差减小时，步长因子的值也会相应减小。

图 10 变步长自适应滤波算法程序设计流程图
采用变步长的自适应滤波算法对二维圆周运动的轨迹进行追踪，滤波结果示于图 11至 13。其中参数。

图 11 X方向自适应滤波结果-变步长自适应滤波

图 12 Y方向自适应滤波结果-变步长自适应滤波

图 13 二维圆周运动轨迹滤波结果-变步长自适应滤波
从X方向与Y方向上的滤波曲线可以看出，变步长的自适应滤波输出结果与期望信号之间的误差更小，固定步长时起始阶段的大幅度波动也消失了，对运动轨迹的追踪效果也更好。
5.解相关自适应滤波
当输入信号之间具有较强的相关性时，自适应滤波的效果并不理想，因此改进自适应滤波算法的一个方法就是消除相邻输入信号序列的相关性，称为解相关自适应滤波。
解相关自适应滤波算法的实现过程为：

该算法通过求解相邻两个输入信号序列的相关系数，在当前输入信号中减去与上一输入信号的相关部分，作为当前的输入信号，实现解相关的自适应滤波。图 14给出了解相关自适应滤波算法的程序设计流程。

图 14 解相关自适应滤波算法流程图
将该算法应用于二维圆周运动的轨迹追踪问题，所得结果示于图 15至 17。

图 15 X方向自适应滤波结果-解相关自适应滤波

图 16 Y方向自适应滤波结果-解相关自适应滤波

图 17 二维圆周运动轨迹滤波结果-解相关自适应滤波
图 15， 16， 17显示了应用解相关自适应滤波算法对二维圆周运动轨迹进行滤波后的结果。X与Y方向上的信号均与期望信号符合的很好，并且最小均方误差的变化曲线也呈现较快的收敛趋势。在二维轨迹图上，滤波轨迹的波动性大大降低，仅在初始阶段存在轻微的波动，但总体上取得了理想的滤波结果，能够满足实际应用的需求。

6.变换域自适应滤波
从解相关自适应滤波算法结果看出，如果能够消除输入信号的相关性，自适应滤波的效果将得到极大的改进，在此基础上，有发展出了变换域的自适应滤波算法。其基本思想是使用一组正交基，将时域信号变换到对应的变换域上，则在变换域上，信号的相关性就会降低，对信号进行归一化后，自相关矩阵特征值的分散度就会降价，从而提高算法的收敛性。基本的变换包括频率域变换，余弦变换，小波变换，分数阶Fourier变换。
（1） 基于频域的自适应滤波
将输入信号和期望信号分别形成N点数据块，然后做N点离散Fourier变换，权系数每N个样点更新一次。对信号进行变换与反变换时，可以利用快速Fourier正变换与逆变换算法，能够有效提高运算速度。
（2） 基于余弦变换域的自适应滤波算法
余弦变换能够较好地近似理想正交变换，基于余弦变换域的LMS自适应滤波算法不仅减小了输入信号的自相关程度，明显提高了收敛速度，减小了权失调噪声，而且该算法的计算量也大大减小。
（3） 基于小波变换域的自适应滤波算法
对自适应滤波器的输入信号进行正交变换，利用小波的时频局部特性，将输入向量正交分解到多尺度空间 。减小了自适应滤波器输入向量自相关阵的谱动态范围，大大增加了算法的收敛步长，提高了收敛速度和稳定性。
（4） 基于分数阶Fourier域的自适应滤波算法
分数阶Fourier变换是一种时频分析工具和旋转算子，信号在分数Fourier域上的表示同时融合了信号在时域和频域的信息。基于分数阶傅里叶变换域的自适应滤波利用前一时刻已获得的滤波器参数等结果，自动调节现时刻的滤波器参数，以适应信号和噪声未知的或随时间变化的统计特性，从而实现最优滤波。

图 18 变换域自适应滤波算法流程图
参考文献
[1]李方伟,张浩. 一种新的变步长LMS自适应滤波算法及其仿真[J]. 重庆邮电大学学报(自然科学版),2009,(05):591-594.
[2]齐林,周丽晓. 变换域自适应滤波算法的研究[J]. 郑州大学学报(理学版),2007,(01):61-66.
[3]冯存前,张永顺. 变步长频域快速自适应收发隔离算法研究[J]. 电子对抗技术,2004,(05):22-25+45.
[4]Deherty J, Porayath R. A robust echo canceler for acoustic environments. IEEE Trans. Circuits and Systems, II, 1997, 44:389-398.
[5]张贤达. 现代信号处理（第三版）.北京:清华大学出版社,2015.
[6]高西全,丁玉美. 数字信号处理-时域离散随机信号处理. 西安:西安电子科技大学出版社, 2002.

代码：
自适应滤波算法：
%该程序实现对二维圆周运动轨迹的自适应滤波
%该程序为主函数，调用不同的子函数实现不同的滤波方法
%子函数1：fun_fplms_filter2-固定步长自适应滤波
%子函数2：fun_cplms_filter2-变步长自适应滤波
%子函数3：fun_lms_filter_der2-解相关自适应滤波
clear
close all
N=2000;
theta=linspace(0,2*pi,N); %极坐标参数
e_x=cos(theta); %x，y方向上的期望信号
e_y=sin(theta);

no_x=normrnd(0,sqrt(0.08),1,N); %高斯白噪声
no_y=normrnd(0,sqrt(0.12),1,N);
m_x=e_x+no_x; %观测信号
m_y=e_y+no_y;
%fixed step
% [Err_x,f_x]=fun_fplms_filter2(e_x,m_x,N,10);
% [Err_y,f_y]=fun_fplms_filter2(e_y,m_y,N,10);
%changed step
% [Err_x,f_x]=fun_cplms_filter2(e_x,m_x,N,10);
% [Err_y,f_y]=fun_cplms_filter2(e_y,m_y,N,10);
%decorrelation
[Err_x,f_x]=fun_lms_filter_der2(e_x,m_x,N,10);
[Err_y,f_y]=fun_lms_filter_der2(e_y,m_y,N,10);

figure
plot(e_x,e_y,'k','linewidth',2)
hold on
plot(m_x,m_y,'b')
hold on
plot(f_x,f_y,'r-')
title('LMS自适应滤波圆周运动轨迹追踪')
legend('期望轨迹','观测轨迹','滤波轨迹')
figure
subplot(211)
plot(e_x,'k')
hold on
plot(m_x,'b')
hold on
plot(f_x,'r')
title('x方向上信号滤波效果对比')
legend('期望信号','观测信号','滤波信号',4)

subplot(212)
plot(Err_x,'k')
title('x方向上滤波方差变化曲线')

figure
subplot(211)
plot(e_y,'k')
hold on
plot(m_y,'b')
hold on
plot(f_y,'r')
title('y方向上信号滤波效果对比')
legend('期望信号','观测信号','滤波信号',4)

subplot(212)
plot(Err_y,'k')
title('y方向上滤波方差变化曲线')

function [SE,x_f]=fun_fplms_filter2(x0,xm,n,m)
%this function conducts the adaptive filtering with fixed step length
x_e=x0;
x_m0=xm;
N=n;
M=m;

x_f=x_m0; %order of filter and initial weight values
w=zeros(1,M);
SE=zeros(1,N);

%% fundmental LMS adptive filter @ Modern SP Zxd P183
rxx=xcorr(x_m0)/N;
Rxx=toeplitz(rxx(N:end));
mui_max=1/max(eig(Rxx));
% trace(Rxx)
% mui_max=1/trace(Rxx); %convergence condition
mui=0.6*mui_max; %initial step length

%% normallized LMS @ Modern SP Zxd P183
% alpha=0.8;beta=2;

%% the iterative filter

x_m=[zeros(1,M) x_m0];
for i=1:N
u_in=x_m(M+i:-1:i+1);
u_out=sum(u_in.*w);
err=x_e(i)-u_out;
% mui=alpha/(beta+sum(u_in.^2));
% mui=0.06;
w=w+mui*u_in*err;
x_f(i)=u_out;
se=x_e-x_f;
SE(i)=sum(se.^2)/N;
end

function [SE,x_f]=fun_cplms_filter2(x0,xm,n,m)
%this function conducts the adaptive filtering with varing length
x_e=x0; %parameter in function mode
x_m0=xm;
N=n;
M=m;

x_f=x_m0; %order of filter and initial weight values
w=zeros(1,M);
SE=zeros(1,N);

%% fundmental LMS adptive filter @ Modern SP Zxd P183
% rxx=xcorr(x_m0)/N;
% Rxx=toeplitz(rxx(N:end));
% mui_max=1/max(eig(Rxx));
% % trace(Rxx)
% % mui_max=1/trace(Rxx); %convergence condition
% mui=0.6*mui_max; %initial step length

%% normallized LMS @ Modern SP Zxd P183
% alpha=0.8;beta=2;
alpha=-6.6;beta=0.18;a=2; %Sigmoid fucntion
%% the iterative filter

x_m=[zeros(1,M) x_m0];
for i=1:N
u_in=x_m(M+i:-1:i+1);
u_out=sum(u_in.*w);
err=x_e(i)-u_out;
mui=beta*(1-exp(alpha*err^a));
% mui=alpha/(beta+sum(u_in.^2));
w=w+mui*u_in*err;
x_f(i)=u_out;
se=x_e-x_f;
SE(i)=sum(se.^2)/N;
end
% x_f(M:M+3)=x_m(M:M+3);% to improve the initial steps of filter
% Err=x_f-x_e;
% Re=Err./x_e;
% find(abs(Re)==max(abs(Re)))
%% plot the filter outputs
figure
plot(x_e,'k')
hold on
plot(x_m0,'b')
hold on
plot(x_f,'r')
% legend('expected','measured','filtered')
legend('期望信号','观测信号','滤波信号')
% figure
% plot(SE,'k')
% title('自适应滤波方差变化曲线')

function [SE,x_f]=fun_lms_filter_der2(x0,xm,n,m)
%this function conducts the adaptive filtering using de-corellation method
N=n;
M=m;
x_e=x0;
x_m0=xm;

x_f=x_m0;
SE=zeros(1,N);
%order of filter and initial weight values
w=zeros(1,M);

%% fundmental LMS adptive filter @ Modern SP Zxd P183
% rxx=xcorr(x_m)/N;
% Rxx=toeplitz(rxx(N:end));
% mui_max=2/trace(Rxx); %convergence condition
% mui=0.4*mui_max; %initial step length

%% normallized LMS @ Modern SP Zxd P183
% alpha=0.5;beta=1;
%% decorrelation LMS @ Modern SP Zxd P183
%to cancle the correlation between u(n) and u(n-1)
rho=0.8;

%% the iterative filter
x_m=[zeros(1,M) x_m0];
u_in0=x_m(M+1:-1:2);
SE(1)=sum((x_e-x_f).^2)/N;
for i=2:N
u_int=x_m(M+i:-1:i+1);
err=x_e(i)-sum(u_int.*w);
corr=sum(u_int.*u_in0)/sum(u_in0.*u_in0); %correlation cofficient
v=u_int-corr*u_in0;
mui=rho*err/sum(u_int.*v);
w=w+mui*v;
u_out=sum(u_int.*w);
x_f(i)=u_out;
u_in0=u_int;
se=x_e-x_f;
SE(i)=sum(se.^2)/N;
end