现代数字信号处理第二章2017

1、1第二章第二章 维纳滤波与卡尔曼滤波维纳滤波与卡尔曼滤波2.1 Wiener filter 介绍2.2 Wiener filter 时域解2.3 Wiener filter Z 域解2.4 Wiener filter 预测2.5 Kalman filter 232.1 Wiener filter 介绍介绍v 滤波与预测滤波与预测 滤波定义滤波定义 所谓滤波所谓滤波,是指在含噪信号是指在含噪信号x(k)=s(k)+v(k) 或其矢或其矢 量信号量信号x(k)=s(k)+v(k)中尽可能排除噪声中尽可能排除噪声v(k)或或v(k) 干扰，而将有用信号干扰，而将有用信号s(k)或或s(k)分离或提取

2、出来。分离或提取出来。 滤波、预测与平滑滤波、预测与平滑 设基于观测过程设基于观测过程x(k)或矢量观测过程或矢量观测过程x(k)，对，对 s(k+)或或s(k+)作最优估计，那么作最优估计，那么 若若=0，就是滤波问题。，就是滤波问题。 若若0，就是预测问题。，就是预测问题。 若若0的约束，非实时： IIR wienerl实时情况则用长为N的有限长度序列来逼近 计算量大！计算量大！增加h(n)的长度来提高精度， FIR wiener11 xsoptxxzHzz 101NmNiiiy ns nh m x nms nh x10,1, 2,0NiijijEsh xxjNEex xs1xxxxxso

3、pthhhopfwienerh2.3 wiener 滤波器滤波器Z域解域解l12 A z w n S z W z s n维纳滤波器输入维纳滤波器输入输出关系输出关系13 A z w n W z v z v n x n s n S z X z a b X z W z x n w n B z利用白化利用白化 的的方法求解维纳方法求解维纳霍霍夫方程夫方程14 H z y ns n x ns nv n( )a 1B z G z x n w n y ns n X z W z( )b x nl15噪声密度越大，噪声密度越大， H越小越小16jwoptHejwoptHe ssPw wPw0wjwoptHe

4、ssPw wPw 0,0,01,0,0,0,0sswjwoptsswsssswsswP wP wH eP wP wP wP wP wP w P wIIR 因果因果wiener 滤波器滤波器17 21和非因果一样，只是多xsoptwzHzB z B z18)8 .01 (5 .01)(B)8 .01)(8 .01 (8 .08 .026 .1)8 .01)(8 .01 ()1)(1 ()(6 .15825 .0)1)(1 (8 .08 .02)8 .01)(8 .01 (8 .08 .021)8 .01)(8 .01 (36.0)()()(2)()()()()( s)(IIR)()(10)(),

5、()()(),(,)8 .0-1)8 .0-136.0)()(111111122121111)(1zzzzzzzzzfzfzzSffzfzzzzzzzzzzSzSzSzSzSzSzSnnxnsnvnvnvnsnxnvzznsnsxxvvsssvvvssvsvsxx得或得令：谱及谱分解求测量数据序列的功率解：的线性最佳估计。得到进行处理维纳滤波器，由它对不相关，试计算一因果与的白噪声，方差为是均值为测量数据加性噪声测量时该信号混入（的功率谱为已知信号19)()(S)(1)(H)8 .01 (6 .0)()(S)5 .01 (3)8 .01 (6 .0)5 .01 ()8 .01 ()8 .01)

6、(8 .01 (36.0)()(S)(S)(S)(S)(S)(S)(/)(S1211111)(1zBzzBzzzBzzzzzzzzzBzzzzzzzBzxsxsxsssvssssvsxsxx传输函数进行因果和逆因果分解对2.4 Wiener 预测预测l惯性： x(n)各点相关，形成惯性 预测可能l一步预测2021211) )()()()()( )()()( )(ypmpmpmpmpmmnxanxEneEmnxanxmnxmhnxn相应误差： 21xsoptwzHzB z B z NxsoptxxzzHzzN步纯预测步纯预测 v(n)=0l最小均方差随N增大而增大lW(n)对X(N)(s(n)的

7、影响就统计平均为0，估算时只需考虑B(z)的惯性，默认有最小均方误差。21 212111NwNoptwzB z B zHzz B zB zB zB z)()(E1022min2nbNneNnw2221211) )()()()()( )()()( )(ypmpmpmpmpmmnxanxEneEmnxanxmnxmhnxn相应误差：一步时域预测：一步时域预测：23自相关，往往已知。只需的互相关，往往未知，和已知估计估计对比：)(xyule)()()(xYule),(wi )(R)()0(R)(1.2 , 1)(R)(-)(R)(R)0(R)(ker.2 , 1)(R)(R1ptssmin21 -N

8、01min21nnxnsWienernnsenermmhneEhopfwienerNlmjmhjmlaneEwalYuleplmlalxspmoxxmoptxxxxpmpmxxxxpmpmxxWiener 判定判定24)(R)()0(R)(1.2 , 1)(R)()(R)(R)0(R)(ker.2 , 1)(R)(Rwiener)(R)()0(R)(1.1 , 0)j (R)(R)(R)()0(R)(0)(R)()(R1ptssmin21 -N01min21x1ptssmin210s1ptssmin21smmhneEhopfwienerNlmjmhjmlaneEwalYuleplmlalmmh

9、neENjmhjmmhneEjmjmhjxspmoxxmoptxxxxpmpmxxxxpmpmxxspmoxxNmoptxxspmoxxmoptx噪声与信号不相关有限脉冲因果2526NiNiiiftxIFNftwftSNftwNftSftxftxIFNftwftSNftwNftSftxftftxftwNftSftSft11iii),(1),(),(COV1),(1),(),(),(1),(),(COV1),(1),(),(),(x),(),(1),(),(),(hwiener修正：和分别对验基于样本平均的时频后27TFR of each HRECG recordsEnsemble avg. o

10、f each TFRCompute ensemble avg.TFR of ensemble avg.Enhance TFRs to attenuate Var & Interference termsCal. H(t,f) to apply TFR of ensemble avg.Reconstruct time domain signal from modiefied linear TFR28292.5 Kalman 滤波器滤波器 l以发明者Rudolf.E.Kalman 而命名，但是在Kanlman之前，Thorvald Nicolai Thiele和Peter Swerling

11、 已经提出了类似的算法lStanley Schmidt 首次实现了Kalman滤波器。在一次对NASA Ames Research Center访问中，卡尔曼发现他的方法对于解决阿波罗计划的轨迹预测很有用，后来阿波罗飞船导航 电脑就使用了这种滤波器。lSwerling(1958),Kalman(1960),Kalman和Bucy(1961)发表的论文。 l 时域滤波器3031问题定义问题定义l系统状态无法直接测量l需要从测量模型中取得最佳估计Measuring DevicesEstimatorMeasurementError SourcesSystem State (desired but n

12、ot known)External ControlsObserved MeasurementsOptimal Estimate of System StateSystemError SourcesSystemBlack Box32l船位置不知l位置函数： y(t)l速度恒定l测量分布符合高斯分布y33010203040506070809010000.020.040.060.040.16T1测量: Mean = z1 and Variance = z1位置最佳估计: (t1) = z1估计方差: 2x (t1) = 2z1船t2 时的位置- Predicted positi

13、on is z1010203040506070809010000.020.040.060.040.16Conceptual Overview-(t2)测量值t2: Mean = z2 and Variance = z2根据测量值修正预测(t2)预测 -(t2)测量z(t2)010203040506070809010000.020.040.060.040.16修正所得即位置最佳估计新方差小于预测及测量方差测量z(t2)修正最佳估计(t2)预测 -(t2)根据之前数据预测- -, - 测量 zk, z最佳估计() = 预测 + (Kalman Gain)

14、 * (测量 预测)方差估计= 方差预测* (1 Kalman Gain)010203040506070809010000.020.040.060.040.16t3, 船移动速度为dy/dt=u(t2)Nave Prediction -(t3)010203040506070809010000.020.040.060.040.16dy/dt = u + w(t2)Nave Prediction -(t3)预测 -(t3)010203040506070809010000.020.040.060.040.16 测量t3修正预测重复前面

15、步骤预测 -(t3)测量 z(t3)修正最佳预测 (t3)40测量仪器测量仪器Estimator测量误差源系统状态外部控制测量观察值Optimal Estimate of System State系统误差源系统系统黑盒子41问题小结问题小结l起始条件(k-1 and k-1)l预测(-k , -k)l用起始条件和模型(例如匀速率)作预测l 测量 (zk)l 修正 (k , k)l用测量值修正预测l最佳估计Kalman 滤波器滤波器42filter)1| 1( )() 1-| 1-( )|( )() 1-| 1-( )|( IIRWIENER) 1-() 1-| 1-( )()|( kalman

16、nnsacnxGnnsannsnGxnns fnnsnsnnsnsnnns可以改为：滤波器的差分方程：则递归实现的因果的最佳线性估计对：时刻及以前所有数据的最佳线性估计时刻及以前所有数据对：应用应用1l假设房间的真实温度为25度，模拟了200个测量值输入，测量值的平均值为25度，但是加入了标准偏差为几度的高斯白噪声（在图中为蓝线）。l为了令卡尔曼滤波器开始工作，设卡尔曼两个零时刻的初始值，是X(0|0)和P(0|0)。因为随着卡尔曼的工作，X会逐渐的收敛。但是对于P，一般不要取0，因为这样可能会令卡尔曼完全相信你给定的X(0|0)是系统最优的，从而使算法不能收敛。设X(0|0)=1度，P(0|

17、0)=10。 l该系统的真实温度为25度，图中用黑线表示。图中红线是卡尔曼滤波器输出的最优化结果（该结果在算法中设置了Q=1e-6，R=1e-1）。43l附matlab下面的kalman滤波程序： lclear lN=200; lw(1)=0; %w为过程噪声 lw=randn(1,N) lx(1)=25; la=1; %a为方程中A(k) lfor k=2:N; lx(k)=a*x(k-1)+w(k-1); lend lV=randn(1,N); %V为观察噪声 lq1=std(V); lRvv=q1.2; lq2=std(x); lRxx=q2.2; lq3=std(w); lRww=q3

18、.2; lc=0.2; %c为方程中H(k) lY=c*x+V; %Y为观察值 lp(1)=0; ls(1)=0; lfor t=2:N; lp1(t)=a.2*p(t-1)+Rww; %p1为方程中p lb(t)=c*p1(t)/(c.2*p1(t)+Rvv); ls(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1); lp(t)=p1(t)-c*b(t)*p1(t); lend lt=1:N; lplot(t,s,r,t,Y,g,t,x,b); 4445residualan heavily th more prediction weightsand increasesK

19、zero todecreasesk -P covarianceerror Predictionpredictionabout sure are weIf,)()()(G)(G) 1-|( )|( ),1-|( G)4) 1|1( )() 1|( )()(a(3) 1|1( ) 1-|( ) 1|( )() 1-|( )2) 1-|1-( ) 1-|( )() 122进行的修正越多测越不准确，利用新息最佳加权系数越大，预预测误差功率越大，修正后得值对新息加权，修正预测选择新息）预测误差预测对测量方程预测cnPRcnPcRncPnannsnnsnnsnnsacnxnnxnxnnnsacnnscnn

20、xnxnnsnnsannsnsnnn460, 122)() 1 | 1 ( ),1 ()0|0( )0()()1 ()()()() 1()( ) 1| 1( )() 1-| 1-( )|( 逐渐下降，直到稳定于可计算和由初始滤波器重要公式总结nsGPsnPcGnnPcRncPGQnanPnnsacnxGnnsannsKalmannn47Blending Factor If we are sure about measurements: Measurement error covariance (R) decreases to zero K decreases and weights resid

21、ual more heavily than prediction If we are sure about prediction Prediction error covariance P-k decreases to zero K increases and weights prediction more heavily than residual48Kalman Kalman 滤波滤波l特点 : 1.数学公式用状态空间描述2.解是递推计算，与wiener滤波不一样，期望响应未知3.采用最小二乘法49-k = Ayk-1 + BukP-k = APk-1AT + Q预测 (Time Upda

22、te)(1) Project the state ahead(2) Project the error covariance ahead修正 (Measurement Update)(1) 计算卡尔曼增益(2) 根据测量zk值更新估计(3) 更新误差方差k = -k + K(zk - H -k )K = P-kHT(HP-kHT + R)-1Pk = (I - KH)P-kl一个是n-1对n时刻估计值，一个是n时刻的测量值，估计值和测量值都存在误差，且误差都假设满足独立的高斯分布lKalman滤波器就是充分结合了估计值和测量值得到n时刻更接近真值的估计结果lKalman滤波器引入状态空间的目的

23、是避免了“像Wiener滤波器一样需要对过去所有0,n-1时刻协方差先验知识都已知”，而直接可以通过上一时刻即n-1时刻的状态信息和均方误差信息就可递推得到n时刻的估计。ln-1对n时刻的估计实际上使用到了所有前0,n-1时刻的信息，信息一直通过最小均方误差进行传递到n-1时刻。基于此，Kalman滤波也需要先验知识，即-1时刻的初始值。50l假设状态空间的n-1时刻估计值和观测空间的n时刻测量值都满足独立高斯分布，Kalman滤波器就是通过高斯分布的乘积运算将估计值和测量值结合，获得最接近真值的n时刻估计。l高斯分布乘积运算的结果仍为高斯分布，高斯分布的均值对应n时刻的估计值，高斯分布的方差

24、对应n时刻的均方误差。51520102030405060708090100-0.7-0.6-0.5-0.4-0.3-0.2-0.10530102030405060708090100-0.7-0.6-0.5-0.4-0.3-0.2-0.10Quick Example Constant ModelLarger value of R the measurement error covariance (indicates poorer quality of measurements)Filter slower to believe measurements slower convergence54Ka

25、lman 滤波本质滤波本质l并不直接估计观测信号， 先估计状态信号l获取wiener解的第归过程l从某个初始状态启动,经迭代运算，最终达到稳态wiener状态55v随机信号的最优预测和滤波随机信号的最优预测和滤波v最优滤波理论与维纳滤波器最优滤波理论与维纳滤波器v横向横向LMSLMS自适应数字滤波器自适应数字滤波器v横向横向RLSRLS自适应数字滤波器自适应数字滤波器v自适应格型滤波器自适应格型滤波器v快速横向滤波（快速横向滤波（FTFFTF）自适应算法）自适应算法v无限脉冲响应自适应滤波器无限脉冲响应自适应滤波器 v盲自适应信号处理盲自适应信号处理v自适应滤波器应用自适应滤波器应用 56自适

26、应信号处理分类：自适应信号处理分类： 自适应天线及自适应天线及 自适应滤波自适应滤波器器 自适应滤波器包括两个过程：滤波过程和自适应过程。自适应滤波器包括两个过程：滤波过程和自适应过程。 此仅考虑后者此仅考虑后者,即滤波器的自适应实现问题；且主要考虑即滤波器的自适应实现问题；且主要考虑 FIR滤波器的自适应实现，其关键是自适应算法。滤波器的自适应实现，其关键是自适应算法。 FIR滤波器的自适应实现指的是：滤波器的自适应实现指的是：M 阶阶FIR滤波器的抽滤波器的抽 头权系数头权系数w1,wM-1可以根据估计误差可以根据估计误差e(n)的大小自动调的大小自动调 节，使得误差在某个统计最优准则下最

27、小。节，使得误差在某个统计最优准则下最小。 滤波器设计最常用的准则是滤波器设计最常用的准则是MMSE准则，即是使滤波器准则，即是使滤波器 实际输出实际输出y(n)与期望响应期望响应d(n)之间的均方误差之间的均方误差 最小最小; 最终达到最终达到Wiener解。解。)(2neE 三、横向LMS自适应数字滤波器v 基本原理基本原理 自适应滤波原理自适应滤波原理l自适应预测可用于语音编码，谱估计， 谱线增强，信号白化57自适应算法参数可调数字滤波器x(n)d(n)y(n)e(n)+-图图3.1 自适应滤波器原理图自适应滤波器原理图自适应建模自适应建模l正向建模正向建模l自适应控制系统，数字滤波器设

28、计，相干估计自适应控制系统，数字滤波器设计，相干估计和地球物理和地球物理l逆向建模逆向建模l自适应控制，语音分析自适应控制，语音分析, 信道均衡，数字滤波器信道均衡，数字滤波器等5859未知系统未知系统自适应处理器x(n)d(n)y(n)+-e(n)+噪声N(n)设备自适应处理器Zy(n)x(n)噪声N(n)d(n)=x(n- )+-e(n)设备图图3.3 自适应建模自适应建模（a）正向建模正向建模 （b）逆向建模逆向建模（a）（b）自适应抵消原理自适应抵消原理60d(n)e(n)+-( )N nx(n)+N(n)( )( )y nN n自适应处理器图图3.4 自适应干扰抵消原理图自适应干扰抵

29、消原理图61e(n)+-12L固定目标信号滤波器多输入自适应处理器+( )( )y nN nx(n)+N(n)传感器阵列延时器图图3.5 自适应阵列信号处理系统原理图自适应阵列信号处理系统原理图多输入干扰抵消。多输入干扰抵消。 E(n) x(n); y(n) N(n)自适应滤波器设计自适应滤波器设计l首先确定滤波器结构-FIR ？ IIR ？格型？l设计自适应算法调整滤波器参数l使得某一特定的代价函数最小62自适应线性组合器自适应线性组合器63自适应处理器0( )x n1( )x n( )Lxn0( )w n1( )w n( )Lw n+y(n)e(n)+-d(n)输入信号矢量l自适应线性组合

30、器： 参数可调的FIR l输入信号X(n)的L+1个元素可以是同一时刻不同信号源所得，也可是同一信号源以前L+1时刻得到l固定权系数lY(n)等于输入矢量x(n)的各元素线性加权l调整权系数的过程叫自适应过程l自适应线性组合器输出就不再是输入信号的线性函数，还包括误差信号e(n)6465min)()()()()()()()()()()()(),(),(),()()()()(22100neEnnxnndnnxndnyndnennnnnnnxnyTTTLLkkk调整权矢量误差信号均值最小准则自适应线性组合器按照字滤波器实际是一个时变横向数多输入情况：wwwwwwww66 式中式中w(n)为第为第n

31、步迭代步迭代(亦即时刻亦即时刻n)的权向量的权向量， 为第为第n步步 迭代的更新步长迭代的更新步长，v(n)为第为第n步迭代的更新方向步迭代的更新方向(向量向量) 三、横向LMS自适应数字滤波器v 基本原理（续）基本原理（续） 最广泛使用的自适应算法是最广泛使用的自适应算法是“下降算法下降算法”)()(2) 1()()() 1()(nxnenwnnnnvww)(n 下降算法的两种实现方式下降算法的两种实现方式 - 自适应梯度算法：自适应梯度算法：LMS算法及其改进算法算法及其改进算法 - 自适应高斯牛顿算法：自适应高斯牛顿算法：RLS算法及其改进算法算法及其改进算法 本节介绍本节介绍LMS类算

32、法，下一节介绍类算法，下一节介绍RLS类算法。类算法。67+1z1z+21( )x n2( )x n0( )w n1( )w n+-e(n)y(n)图图3.13 自适应线性组合器采用自适应线性组合器采用LMS算法的计算流程图算法的计算流程图lLMS算法核心思想： 用平方误差代替均方误差l权系数调整有噪声-权矢量在最近权矢量附近随机起伏。l收敛：每次迭代调整权矢量，获得多个x(n),并统计平均根据梯度平均值来调整权矢量，则必然得到理想的最佳权矢量6869表示滤波器在表示滤波器在n时刻的估计误差，并定义均方误差时刻的估计误差，并定义均方误差: 三、横向LMS自适应数字滤波器v 基本基本LMSLMS

33、算法算法算法推导算法推导 )()()(nndneHxw22)()(E21)(E21)(nndnenJHxw为代价函数；它相对于滤波器抽头权向量为代价函数；它相对于滤波器抽头权向量w的梯度为的梯度为 )()()(E)()(E)()()()(E)(*nnnndnnnndnnJHHwxxxwxx它是式它是式(4)的向量形式。若上式中的数学期望项用它的的向量形式。若上式中的数学期望项用它的瞬时值代替瞬时值代替, 即得真实梯度向量的估值即得真实梯度向量的估值(瞬时梯度瞬时梯度 )：)13()()()()()()()()()()()()(*nennnndnnnnndnnJHHxwxxwxxx令令70其中其

34、中 三、横向LMS自适应数字滤波器v 基本基本LMSLMS算法算法算法推导算法推导(续续) 设用最速下降法更新滤波器权向量，则有如下算法：设用最速下降法更新滤波器权向量，则有如下算法：)14()()()()()()()() 1(*annennnJnnnxwww)14()(*)()()()()()(bnnndnnndneTHwxxw式式 (14)所示算法就是著名的最小均方自适应算法，简所示算法就是著名的最小均方自适应算法，简称称LMS算法算法。它由。它由Widrow在在20世纪世纪60年代初提出的。年代初提出的。 71 而且可以证明而且可以证明LMS自适应滤波器的权向量收敛于维纳解：自适应滤波器的权向量收敛于维纳解：三、横向LMS自适应数字滤波器v 基本基本LMSLMS算法算法R算法收敛性算法收敛性 容易验证，瞬时梯度向量是真实梯度向量的无偏估计：容易验证，瞬时梯度向量是真实梯度向量的无偏估计： )()()()()(E)(E*nJnnndnnJHwxxrRw1)1(limkEk 条件是条件是 max/20 LMS算法还