维纳滤波和卡尔曼滤

1、 随机信号或随机过程随机信号或随机过程(random process)是普遍存在的。是普遍存在的。 通常把对信号或系统功能起干扰作用的通常把对信号或系统功能起干扰作用的随机信号称之为噪声。随机信号称之为噪声。 噪声按功率谱密度划分可以分为白噪声噪声按功率谱密度划分可以分为白噪声（white noise）和色噪声（）和色噪声（color noise） 均值为均值为0的白噪声叫纯随机信号（的白噪声叫纯随机信号（pure random signal）。）。 任何随机信号都可看成是纯随机信任何随机信号都可看成是纯随机信号与确定性信号并存的混合随机信号与确定性信号并存的混合随机信号或简称为随机信号。号或

2、简称为随机信号。 要区别干扰（要区别干扰（interference）和噪声）和噪声( noise)两种事实和两个概念。非目两种事实和两个概念。非目标信号（标信号（nonobjective signal）都）都可叫干扰。可叫干扰。 干扰可以是确定信号，如国内的50Hz工频干扰。干扰也可以是噪声，纯随机信号（白噪声）加上一个直流成分（确定性信号），就成了最简单的混合随机信号。 医学数字信号处理的目的是要提取包含在随机信号中的确定成分，并探求它与生理、病理过程的关系，为医学决策提供一定的依据 例如从自发脑电中提取诱发脑电信号，就是把自发脑电看成是干扰信号，从中提取出需要的信息成分。因此我们需要寻找一

3、种最佳线性滤波器，当信号和干扰以及随机噪声同时输入该滤波器时，在输出端能将信号尽可能精确地表现出来。匹配滤波 匹配滤波着眼点不是尽可能保持信号不是真，而是提高输出的信噪比 维纳滤波和卡尔曼滤波就是用来解决这样一类问题的方法：从噪声中提取出有用的信号。 这种线性滤波方法也被看成是一种估计问题或者线性预测问题。 维纳滤波-对真实信号的最小均方误差(MMSE, minimum mean-square error)估计问题.设有一个线性系统，它的单位脉冲响应是h(n)，当输入一个观测到的随机信号x(n)，简称观测值，且该信号包含噪声w(n)和有用信号s(n)，简称信号，也即 )()()(nwnsnx

4、(9-1)则输出)(ny为mmnxmhnhnxny)()()()()( (9-2)()()(nwnsnx)( )(nsny )(nh图9-1 维纳滤波器的输入输出关系平滑滤波预测 niixinhns0 11npniixinhns 10Niixinhns这里我们主要考虑滤波问题，即线性估计根据其取值范围不同通常有下面几种情况:从图9-1的系统框图中估计到的信号和我们期望得到的有用信号可能不完全相同，这里用来表示真值和估计值之间的误差)( )()(nsnsne (9-3)显然是随机变量，维纳滤波和卡尔曼滤波的误差准则就是最小均方误差准则：22)( )()(nsnsEneE(9-4) 维纳滤波和卡尔

5、曼滤波都是解决线性滤波和预测问题的方法，并且都是以均方误差最小为准则的，在平稳条件下两者的稳态结果是一致的。但是它们解决问题的方法有很大区别。 维纳滤波是根据全部过去观测值和当前观测值来估计信号的当前值，因此它的解形式是系统的传递函数或单位脉冲响应； 卡尔曼滤波是用当前一个估计值和最近一个观测值来估计信号的当前值，它的解形式是状态变量值。 维纳滤波只适用于平稳随机过程，卡尔曼滤波就没有这个限制。 设计维纳滤波器要求已知信号与噪声的相关函数，设计卡尔曼滤波器要求已知状态方程和量测方程，当然两者之间也有联系。第一节第一节 维纳滤波器的时域解维纳滤波器的时域解（Time domain solutio

6、n of the Wiener filter） 设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳霍夫（WienerHopf）方程。我们从时域入手求最小均方误差下的 ， 用 表示最佳线性滤波器。这里只讨论因果可实现滤波器的设计。)(nhopt0n当0,h(n)0mm)h(m)x(n(n)s y(n)(9-5)202) )()()()(mmnxmhnsEneE(9-6)2 , 1 , 00)() )()()(20jjnxmnxmhnsEmopt(9-7)即0)()()()()(0jjnxmnxEmhjnxnsEmopt(9-8)用相关函数R来表达

7、上式，则得到维纳霍夫方程的离散形式：0)()()(0jmjRmhjRmxxoptxs从维纳霍夫方程中解出的h就是最小均方误差下的最佳h， 即)(nhopt 求到 )(nhopt，这时的均方误差为最小： 20min2) )()()()(moptmnxmhnsEneE )()()()()()()(2)(0002mroptoptmrnxrhmnxmhmnxmhnsnsE000)()()()()(2)0(mrxxoptoptmxsoptssrmRrhmhmRmhR由式（9-9）进一步化简得：0min2)()()0()(mxsoptssmRmhRneE(9-10)二、有限脉冲响应法求解维纳霍夫方程二、有

8、限脉冲响应法求解维纳霍夫方程 如何去求解维纳霍夫方程，即式（9-9）中解的问题，设是一个因果序列且可以用有限长（N点长）的序列去逼近它，则式(9-5) (9-10)分别发生变化10)()()( )(Nmmnxmhnsny (9-11)2102) )()()()(NmmnxmhnsEneE (9-12) 12 , 1 , 00)() )()()(210NjjnxmnxmhnsENmopt(9-13)1, 1 , 0)()()()()(10NjjnxmnxEmhjnxnsENmopt .(9-14)1, 2 , 1 , 0)()()(10NjmjRmhjRNmxxoptxs (9-15)于是得到N

9、个线性方程：)0() 1()2() 1 () 1()0() 1(1)2() 1()0() 1 () 1 ()0() 1 (1) 1() 1() 1 () 1 ()0()0()0(0 xxxxxxxsxxxxxxxsxxxxxxxsRNhNRhNRhNRNjNRNhRhRhRjNRNhRhRhRj写成矩阵形式有： ) 1() 1 ()0() 1() 1 ()0()0()2() 1()2()0() 1 () 1() 1 ()0(NRRRNhhhRNRNRNRRRNRRRxsxsxsxxxxxxxxxxxxxxxxxx(9-16)简化形式：RxxH=Rxs (9-17)式中，Hh(0) h(1) h

10、(N-1)，是待求的单位脉冲响应；Rxs) 1(),1 (),0(NRRRxsxsxs，是互相关序列； Rxx )0()2() 1()2()0() 1 () 1() 1 ()0(xxxxxxxxxxxxxxxxxxRNRNRNRRRNRRR，是自相关矩阵 只要Rxx是非奇异的，就可以求到H：H=Rxx1 Rxs (9-18)求得)(nhopt后，这时的均方误差为最小： 210min2) )()()()(NmoptmnxmhnsEneE )()()()()()()(2)(1010102NmNroptoptNmrnxrhmnxmhmnxmhnsnsE101010)()()()()(2)0(NmNr

11、xxoptoptNmxsoptssrmRrhmhmRmhR 10min2)()()0()(NmxsoptssmRmhRneE (9-19)用有限长的)(nh来实现维纳滤波时，当已知观测值的自相关和观测值与信号的互相关时就可以按照式（9-15）在时域里求解 )(nhopt但是当N比较大时，计算量很大，并 且涉及到求自相关矩阵的逆矩阵问题。 )(ns)(nw0)()(mRmRwssw则有)()()()()()()()(mRmnsnwmnsnsEmnsnxEmRssxs)()()()()()()(mRmRmnwmnsnwnsEmRwwssxx则式（9-15）和式（9-19）化为：1, 2 , 1 ,

12、 0)()()()(10NjmjRmjRmhjRwwNmssoptss(9-20) 10min2)()()0()(NmssoptssmRmhRneE (9-21)【例9-1】已知图9-1中 )()()(nwnsnx且 )(ns与 )(nw统计独立，其中 )(ns的自相关序列为 mssmR6 . 0)()(nw， 是方差为1的单位白噪声，试设计一个N=2维纳滤波器来估计 )(ns，并求最小均方误差。 mssmR6 . 0)()()(mmRww，代入式（9-20）得 ) 1 (2)0(6 . 06 . 01) 1 (6 . 0)0(210hhjhhj解得： )0(h0.451， ) 1 (h0.1

13、65。 将上述结果代入式（9-21），求得最小均方误差：45. 0) 1 (6 . 0)0(1)()()0()(10min2hhmRmhRneEmssss若要进一步减小误差可以适当增加维纳滤波的阶数，但相应的计算量也会增加。 解依题意，已知信号的自相关和噪声的自相关为： 三、预白化法求解维纳霍夫方程 从上面分析知求解维纳霍夫方程比较复杂，本节用波德（Bode）和香农（Shannon）提出的白化的方法求解维纳霍夫方程，得到系统函数 H(z)。 随机信号都可以看成是由一白色噪声 w1(n) 激励一个物 理可实现 的系统 或模型的响应，如图9-2所示， 其中A(z) 表示系统的传递函数。由于x(n)

14、 = s (n) + w(n)，在图9-2的基础上给出x(n)的信号模型，图9-3所示。把这两个模型合并最后得到维纳滤波器的信号模型，图9-4所示，其中传递函数用B(z)表示。 图9-2 的信号模型)(ns)(zA)(1nw 图9-3 )(nx的信号模型 )(ns)(ns)(zA)(1nw)(nw)(nx图9-4 维纳滤波器的输入信号模型)(nx)(zB)(1nw白噪声的自相关函数为 ，它的z变换就等于 。图9-2中输出信号的自相关函数为 ，根据卷积性质有)()(2111mmRwww21w)(mRss )()()()()()()(11rkssrmnwraknwkaEmnsnsEmR)()()(

15、11rkmRrakawwkr令 krl上式 lkwwwwkllkakalmRlmRlkaka)()()()()()(1111令 )()()()()(lalalkakalfl)()()()()()()()(111111mamamRmfmRlflmRmRwwwwlwwss对式（9-22）进行Z变换得到系统函数和相关函数的z变换之间的关系：)()()(121zAzAzRwss (9-23)同样，对图9-4进行z变换得)()()(121zBzBzRwxx(9-24)图9-4中利用卷积性质还可以找到互相关函数之间的关系：)()()()()()(1mnsknwkbEmnsnxEmRkxs)()()()(1

16、1mbmRkmRkbswksw两边z变换得到)()()(11zBzRzRswxs (9-25)如果已知观测信号的自相关函数，求它的z变换， 然后找到该函数的成对零点、极点，取其中在单位圆内的那一半零点 、极点构成)(zB，另外在)(zB)(1zB单位圆外的零、极点构成，这样就保证了 是因果的，并且是最小相位系统。从图9-4可得)()(1)(1zXzBzW (9-26)由于系统函数)(zB的零点和极点都在单位圆内，)(1zB)(nx即是一个物理可实现的最小相位系统，则也是一个物理可实现的最小相移网络函数。我们就可以利用式（9-26）对 )(nx)(1nw)(1zB进行白化，即把当作输入，当作输出

17、，是系统传递函数。将图9-1重新给出，待求的问题就是最小均方误差下的最佳 ，如图9-5(a)所示，为了便于求这个 ，将图9-5(a)的滤波器分解成两个级联的滤波器： 和G（z）， 如图9-5(b)所示，则)(zH)(zHopt)(1zB)()()(zBzGzH (9-27)()()(nwnsnx)( )(nsny)(nh（a）)(nx)( )(nsny)(1zB)(zG)(1nw（b）图9-5 利用白化方法求解模型有了上述的模型后，白化法求解维纳霍夫方程步骤如下： 1）对观测信号 的自相关函数 求z变换得到 ； 2）利用等式 ，找到最小相位系统 ； 3）利用均方误差最小原则求解因果的G（z）；

18、 4） ，即得到维纳霍夫方程的系统函数解。)(nx)(mRxx)(zRxx)()()(121zBzBzRwxx)(zB)()()(zBzGzH在上述步骤中， 可以通过已知的观测信号的自相关函数来求得，因而求解 的问题就归结为求解G（z）的问题了。由于G（z）的激励源是白噪声，求解变得容易多了，下面我们分析步骤3的求解过程。 按图9-5(b)有:)(zB)(zHopt01)()()( )(mmnwmgnsny (9-28)均方误差为：2012) )()()()(mmnwmgnsEneE )()()()()()()(2)(0011012mrmrnwrgmnwmgmnwmgnsnsE000)()()

19、()()(2)0(111mrwwmswssrmRrgmgmRmgR由于 ，代入上式，并且进行配方得：)()(2111mmRwww02202)()()(2)0()(11mwmswssmgmRmgRneE02202)(1)()()0(11111mswwmwswwssmRmRmgR均方误差最小也就是上式的中间一项最小，所以0,)()(211mmRmgwswopt (9-30)注意，这里的)(mg是因果的。对该式求z变换，得到211)()(wswoptzRzG (9-31)(1zRsw)(1mRsw表示对求单边z变换。所以维纳霍夫方程的系统函数解表示为：)()()(zBzGzHoptopt)()(21

20、1zBzRwsw由式（9-25）上式可以表示为：)()()(zBzGzHoptopt)()(211zBzRwsw)()(/ )(211zBzBzRwxs因果的维纳滤波器的最小均方误差为： min2)(neE022)(1)0(11mswwssmRRmswwssmumRR)()(1)0(2211 (9-33)min2)(neEcxsoptsszdzzRzHzRj)()()(211 (9-34)()()(nwnsnx)(ns)(nw)(nsmssmR8 . 0)()(nw)(ns【例9-2】已知图9-1中，且统计独立，其中的自相关序列为，是方差为1的单位白噪声，试，并求最小均方误差。与设计一个物理可

21、实现的维纳滤波器来估计mssmR8 . 0)()()(mmRww0)(mRsw)()(mRmRssxs解依题意，已知，步骤1)()()(mRmRmRwwssxx求z变换25. 18 . 0 ,)8 . 01)(8 . 01 ()5 . 01)(5 . 01 (6 . 11)8 . 01)(8 . 01 (36. 0)(111zzzzzzzzRxx步骤2由于)()()(121zBzBzRwxx，容易找到最小相位系统和白噪声方差zzzzB8 . 0 ,8 . 015 . 01)(1125. 1,8 . 015 . 01)(1zzzzB6 . 121w步骤3利用式（932） )(zHopt)()(/

22、 )(211zBzBzRwxs)5 . 01)(8 . 01 (36. 0)5 . 01 (6 . 18 . 01111zzzz对括号里面求反变换，注意括号内的收敛域为28 . 0 z)5 . 01)(8 . 01 (36. 011zzZ) 1()2(6 . 0)()8 . 0(6 . 0nununn取因果部分，也就是第一项，所以118 . 0116 . 0)5 . 01)(8 . 01 (36. 0zzz)(zHopt0,)5 . 0(375. 0)(nnhn步骤4最小均方误差为 min2)(neEcxsoptsszdzzRzHzRj)()()(211cdzzzzzj)5 . 0)(25.

23、1)(8 . 0()5 . 0625. 0(45. 021取单位圆为积分围线，有两个单位圆内的极点，0.8和0.5，求它们的留数和，所以min2)(neE375. 0)25. 15 . 0)(8 . 05 . 0()5 . 05 . 0*625. 0(45. 0)5 . 08 . 0)(25. 18 . 0()5 . 08 . 0*625. 0(45. 0第二节 维纳预测器（Wieners Predictor） 上节讨论的维纳滤波器是一种估计器，是用观测到的当前 和全部过去的数据 、 、来估计当前的 信号值。本节将讨论维纳预测器，它同样也是一种估计器，是用过去的观测值来估计当前的或将来的信号

24、，N ，也是用真值和估计值的均方误差最小为估计准则。)(nx) 1( nx)2( nx)( )(nsny)( )(Nnsny0一、因果的维纳预测器一、因果的维纳预测器 图9-6就是维纳预测器的模型，N0， 是希望得到的输出，而 表示实际的估计值。)(nyd)(ny)()()(nwnsnx)( )(Nnsny)(nh)()(Nnsnyd图9-6维纳预测器本节和上节一样着重讨论预测器的系统函数以及预测的均方误差，维纳预测器和维纳滤波器比较类似，因而分析方法也都可以借鉴前面的内容。 )(nh0, 0)(nnh当对于图9-6模型，设是物理可实现的，也即，则有是因果序列：0)()()( )(mmnxmh

25、Nnsny (9-35)202) )()()()(mmnxmhNnsEneE (9-36)要使得均方误差最小，则将上式对各，m0，1，求偏导，并且等于零，得：)(mh2 , 1 , 00)() )()()(20jjnxmnxmhNnsEmopt (9-37)即 用相关函数R来表达上式：0)()()()()(0jjnxmnxEmhjnxNnsEmopt(9-38)0)()()(0jmjRmhjNRmxxoptxs (9-39)()(Nnsnyd)()()()(NmRNmnsnxEmRxsxyd由于，则，z变换得)()(zRzzRxsNxyd)()(11zRzzRxsNxyd (9-40)()(N

26、nsnyd)()(nsnyd借鉴维纳滤波器的结果类似给出维纳预测器的最佳传递函数，对应维纳预测器，对应维纳滤波器，故因果的预测器的传递函数为： )(zHopt)()(/ )(211zBzBzRwxyd)()(/ )(211zBzBzRzwxsN (9-41)最小均方误差为 min2)(NneE022)(1)0(11mywwssmRRd (9-42)利用帕塞伐尔(Parseval)定理，上式可用z域来表示min2)(NneEcxyoptsszdzzRzHzRjd)()()(211cxsNoptsszdzzRzzHzRj)()()(211)()()(nwnsnx)(ns)(nw)(nsmssmR8

27、 . 0)()(nw) 1( ns【例9-3】已知图7-6中，且与统计独立，其中的自相关序列为，是方差为1的单位白噪声，并求最小均方误差。试设计一个物理可实现的维纳预测器估计mssmR8 . 0)()()(mmRww0)(mRsw解依题意已知， )()(mRmRssxs)()()(mRmRmRwwssxx求z变换：25. 18 . 0 ,)8 . 01)(8 . 01 ()5 . 01)(5 . 01 (6 . 11)8 . 01)(8 . 01 (36. 0)(111zzzzzzzzRxx由于)()()(121zBzBzRwxx，容易找到最小相位系统和白噪声方差：zzzzB8 . 0 ,8

28、. 015 . 01)(1125. 1,8 . 015 . 01)(1zzzzB6 . 121w由式（9-41），N1，)(zHopt)()(/ )(211zBzBzzRwxs)()(/ )(211zBzBzzRwss求z变换：25. 18 . 0 ,)8 . 01)(8 . 01 ()5 . 01)(5 . 01 (6 . 11)8 . 01)(8 . 01 (36. 0)(111zzzzzzzzRxx由于)()()(121zBzBzRwxx，容易找到最小相位系统和白噪声方差：zzzzB8 . 0 ,8 . 015 . 01)(1125. 1,8 . 015 . 01)(1zzzzB6 .

29、121w由式（9-41），N1，)(zHopt)()(/ )(211zBzBzzRwxs)()(/ )(211zBzBzzRwss)5 .01)(8 .01 (36.0)5 .01 (6 .18 .01111zzzzz对括号里面求z反变换，注意括号内的收敛域为：28 . 0 z，)5 . 01)(8 . 01 (36. 011zzzZ) 1()2(2 . 1)()8 . 0(48. 0nununn取因果部分，也就是第一项，所以118.01148.0)5.01)(8.01(36.0zzzz)(zHopt11115 . 013 . 0)8 . 01 (48. 0)5 . 01 (6 . 18 .

30、01zzzz0,)5 . 0(3 . 0)(nnhn把上式写成差分方程形式有：)( 5 . 0)(3 . 0) 1( nsnxnsmin2)(neEcssoptsszdzzRzzHzRj)()()(2111最小均方误差为：cdzzzj6 . 0)8 . 01)(5 . 0(36. 021二、纯预测器（二、纯预测器（N步）步） 纯预测器指的是 0的情况下，对 的预测。如图7-7所示。)(nw)(Nns)()(nsnx)( )(Nnxny)(nh)()(Nnxnyd图9-7 N步纯预测器)()(nsnx)()()(121zBzBzRwxx)()()()()(121zBzBzRzRzRwxsssxx

31、这时，用白化法来求解预测器的系统函数。因为，从而有： (9-44)(zHopt)()(/ )(211zBzBzRzwxsN)()()()(/ )()(211211zBzBzzBzBzBzBzNwwN将上式代入式（9-41）、（9-43）得：min2)(NneEcxsNoptsszdzzRzzHzRj)()()(211cNNwzdzzBzzBzzBzBj)()()()(21121假设B(z）是b（n）的z变换，且b（n）是实序列，则上式可以利用帕塞伐尔定(Parseval)理进一步化简：min2)(NneE )()()()(221nnwNnbnuNnbnb (9-46)又因为B(z）是最小相位系

32、统，一定是因果的，上式可以简化min2)(NneE102202022)( )()(11NnwnnwnbNnbnb (9-47)上式说明最小均方误差随着N的增加而增加，也即预测距离越远误差越大。 )()(nsnx)(nsmssmR8 . 0)()( NnsmssmR8 . 0)(【例9-4】已知图7-7中，其中的自相关序列为，试设计一个物理可，并求最小均方误差。，则实现的维纳预测器来估计解依题意，已知25. 18 . 0 ,)8 . 01)(8 . 01 (36. 0)(1zzzzRss因为 )()()()()(121zBzBzRzRzRwxsssxx容易找到最小相位系统和白噪声方差:zzzB8

33、 . 0 ,8 . 011)(1)(8 . 0)(nunbn25. 1,8 . 011)(1zzzB36. 021w利用式（9-45）:)(zHopt8 . 011)8 . 01 ()()(11zzzzBzBzNN)(8 . 08 . 01111NnuzzZNnN0n因为 ，只取的部分，有:)(8 . 08 . 01111nuzzZNnN)8 . 0118 . 08 . 01111zzzNN)(zHoptNoptzH8 . 0)(回到z域有：，代入得:最小均方误差为：min2)(NneENNnnNnwnunb210210228 . 01)(8 . 036. 0)(1它说明当N越大，误差越大，当

34、N0时，没有误差。把上述结果用模型表示如图9-8所示。)()(nsnx)(8 . 0)( )(nsNnsnyN)(8 . 0)(nnhNopt图9-8 例题9-3的纯预测模型三、一步线性预测器三、一步线性预测器 对于纯预测问题，有 ，然而预测的问题常常是要求在过去的p个观测值的基础上来预测当前值，也就是0)()()( )(mmnxmhNnxnypmmnxmhnxny1)()()( )(这就是一步线性预测公式，常常用下列符合表示pmpmmnxanx1)()( (9-48)式中p为阶数，)(mhapm。预测的均方误差为:) )()()(212pmpmmnxanxEneEpmplxxplpmpmxx

35、pmxxmlRaamRaR111)( )( 2)0( (9-49)要使得均方误差最小，将上式右边对pma求偏导并且等于零，得到p个等式 pmaplmlRalRpmxxpmxx, 2 , 10)()(1 (9-50)最小均方误差：pmxxpmxxmRaRneE1min2)()0()( (9-51)式（9-50）就是YuleWalker（Y-W）方程，把YuleWalker（Y-W）方程和维纳霍夫方程进行比较，维纳霍夫方程要估计的量是s(n)，Y-W方程要估计的量是x(n)本身，因而解维纳霍夫方程要已知x(n)、y(n)的互相关函数，实际中这个互相关函数往往是未知的，而解Y-W方程只需要知道观测信

36、号的自相关函数。因此Y-W方程比W-H方程更具有实用价值。)( nsmssmR8 . 0)(mssxxmRmR8 . 0)()(64. 0)2(, 8 . 0) 1(, 1)0(xxxxxxRRR例9-5】已知图9-7中x(n)=s(n)，其中的自相关序列为的可实现的一步线性预测器，并求最小均方误差。解， ，试设计一个p2利用Y-W方程2 , 10)()(21lmlRalRmxxpmxx，可以列出2个方程式08 . 064. 008 . 08 . 02121ppppaaaa0, 8 . 021ppaa) 1(8 . 0)( nxnx解得：，也即36. 064. 01)()0()(1min2pm

37、xxpmxxmRaRneE结果和例（9-4）N=1时一致。),(),(),(ftWftSftXii),(),(),(1),(1ftWftSftXNftXNiiiNitwNtstxiNi, 2 , 1),(1)()(1对每次观测用短时傅立叶变换求时频表示（TFR）：对N次观测的时频表示（TFR）求平均：，样本平均的时频表示（TFR）为： (1)样本平均为：),(1),(),(ftWNftSftX (2)从式（2）可以得到一个基于样本平均的简单时频平面后验维纳滤波器：),(1),(),(),(ftWNftSftSfth (3) ),(1),(),(1),(),(),(11NiiiNiiiftXIF

38、NftWftSCOVNftWftSftX .(4),(),(),(),(),(),(ftXIFftWftSCOVftWftSftX (5),(ftW),(ftWi式中COV表示信号和噪声之间的方差，也就是考虑了信号和噪声并非相互独立；IF是干扰项；表示样本平均的噪声功率；表示样本噪声功率的平均。 原信号是线性调频信号，观测信号混有白噪声 在图9-10中每一个图中从上至下分别表示：测量的单个样本，样本平均，TFPW滤波器估计的信号，原始信号。图9-10的初始信噪比设为12dB，TFPW与叠加平均法相比，信噪比有5个dB左右的改善。五、需要进一步研究的问题五、需要进一步研究的问题 FPW滤波中由于

39、有二次TFR中的相关噪声以及IF项，滤波器可能包含虚部，也就是包含信号的相位信息，直接在时频平面上考虑相位问题还需要进一步研究。 简单来说，卡尔曼滤波器是一个“optimal recursive data processing algorithm（最优化自回归数据处理算法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。 他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。 近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。引入事列 假设我们要研究的对象是一个房间的温度。 根据你的经验判断，这个房间的温

40、度是恒定的，也就是下一分钟的温度等于现在这一分钟的温度（假设我们用一分钟来做时间单位）。 假设你对你的经验不是100%的相信，可能会有上下偏差几度。 我们把这些偏差看成是高斯白噪声（White Gaussian Noise），也就是这些偏差跟前后时间是没有关系的而且符合高斯分配（Gaussian Distribution）。 另外，我们在房间里放一个温度计，但是这个温度计也不准确的，测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。 对于某一分钟我们有两个有关于该房间的温度值：你根据经验的预测值（系统的预测值）和温度计的值（测量值） 要估算k时刻的是实际温度值。 首先你要根据k-1时刻的

41、温度值，来预测k时刻的温度。因为你相信温度是恒定的，所以你会得到k时刻的温度预测值是跟k-1时刻一样的，假设是23度，同时该值的高斯噪声的偏差是5度（5是这样得到的：如果k-1时刻估算出的最优温度值的偏差是3，你对自己预测的不确定度是4度，他们平方相加再开方，就是5）。 然后，你从温度计那里得到了k时刻的温度值，假设是25度，同时该值的偏差是4度。 用于估算k时刻的实际温度有两个温度值，分别是23度和25度。究竟实际温度是多少呢？相信自己还是相信温度计呢？究竟相信谁多一点，我们可以用他们的covariance来判断。因为Kg2=52/(52+42)，所以Kg=0.78，我们可以估算出k时刻的实

42、际温度值是：23+0.78*(25-23)=24.56度。可以看出，因为温度计的covariance比较小（比较相信温度计），所以估算出的最优温度值偏向温度计的值。 已经得到k时刻的最优温度值了，下一步就是要进入k+1时刻，进行新的最优估算。到现在为止，好像还没看到什么自回归的东西出现。对了，在进入k+1时刻之前，我们还要算出k时刻那个最优值（24.56度）的偏差。算法如下：(1-Kg)*52)0.5=2.35。 这里的5就是上面的k时刻你预测的那个23度温度值的偏差，得出的2.35就是进入k+1时刻以后k时刻估算出的最优温度值的偏差（对应于上面的3）。 就是这样，卡尔曼滤波器就不断的把cov

43、ariance递归，从而估算出最优的温度值。他运行的很快，而且它只保留了上一时刻的covariance。上面的Kg，就是卡尔曼增益（Kalman Gain）。他可以随不同的时刻而改变他自己的值 卡尔曼滤波的算法The Kalman Filter Algorithm1、预估计 X(k|k-1)=A X(k-1|k-1)+B U(k) . (1) 2、计算预估计协方差矩阵 Z(k)=H X(k)+V(k)（2） 3、计算卡尔曼增益矩阵 Kg(k)= P(k|k-1) H / (H P(k|k-1) H + R) . (3) 4、更新估计 X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H

44、 X(k|k-1) (4) 5、计算更新后估计协防差矩阵 P(k|k)=（I-Kg(k) H）P(k|k-1) . (5) 卡尔曼滤波的算法The Kalman Filter Algorithm 首先引入一个离散控制过程的系统。该系统可用一个线性随机差分方程（Linear Stochastic Difference equation）来描述： X(k)=A X(k-1)+B U(k)+W(k) 再加上系统的测量值：Z(k)=H X(k)+V(k) 上两式子中，X(k)是k时刻的系统状态，U(k)是k时刻对系统的控制量。A和B是系统参数，对于多模型系统，他们为矩阵。Z(k)是k时刻的测量值，H是

45、测量系统的参数，对于多测量系统，H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声(White Gaussian Noise)，他们的covariance 分别是Q，R（这里我们假设他们不随系统状态变化而变化）。 首先我们要利用系统的过程模型，来预测下一状态的系统。假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态：X(k|k-1)=A X(k-1|k-1)+B U(k) . (1)式(1)中，X(k|k-1)是利用上一状态预测的结果，X(k-1|k-1)是上一状态最优的结果，U(k)为现在状态的控制量，如果没有控制量，它可以为0 到现在为

46、止，我们的系统结果已经更新了，可是，对应于X(k|k-1)的covariance还没更新。我们用P表示covariance：P(k|k-1)=A P(k-1|k-1) A+Q (2)式(2)中，P(k|k-1)是X(k|k-1)对应的covariance，P(k-1|k-1)是X(k-1|k-1)对应的covariance，A表示A的转置矩阵，Q是系统过程的covariance。Q是系统过程的covariance。式子1，2就是卡尔曼滤波器5个公式当中的前两个，也就是对系统的预测。 Kg为卡尔曼增益(Kalman Gain)： Kg(k)= P(k|k-1) H / (H P(k|k-1) H

47、 + R) (3) 有了现在状态的预测结果，然后再收集现在状态的测量值。结合预测值和测量值，可以得到现在状态(k)的最优化估算值X(k|k)： X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H X(k|k-1) .(4)其中 到现在为止，已经得到了k状态下最优的估算值X(k|k)。但是为了要令卡尔曼滤波器不断的运行下去直到系统过程结束，还要更新k状态下X(k|k)的covariance：P(k|k)=（I-Kg(k) H）P(k|k-1) (5) 简单例子（A Simple Example） 把房间看成一个系统，然后对这个系统建模。当然，我们建的模型不需要非常地精确。我们所知道的这

48、个房间的温度是跟前一时刻的温度相同的，所以A=1。没有控制量，所以U(k)=0。因此得出：X(k|k-1)=X(k-1|k-1) . (6)式子（2）可以改成：P(k|k-1)=P(k-1|k-1) +Q (7) 因为测量的值是温度计的，跟温度直接对应，所以H=1。式子3，4，5可以改成以下： Kg(k)= P(k|k-1) / (P(k|k-1) + R) (8) X(k|k)= X(k|k-1)+Kg(k) (Z(k)-X(k|k-1)(9)P(k|k)=（1-Kg(k)）P(k|k-1) .(10) 现在我们模拟一组测量值作为输入。假设房间的真实温度为25度，我模拟了200个测量值，这些

49、测量值的平均值为25度，但是加入了标准偏差为几度的高斯白噪声（在图中为蓝线）。 为了令卡尔曼滤波器开始工作，我们需要告诉卡尔曼两个零时刻的初始值，是X(0|0)和P(0|0)。他们的值不用太在意，随便给一个就可以了，因为随着卡尔曼的工作，X会逐渐的收敛。但是对于P，一般不要取0，因为这样可能会令卡尔曼完全相信你给定的X(0|0)是系统最优的，从而使算法不能收敛。我选了X(0|0)=1度，P(0|0)=10。该系统的真实温度为25度，图中用黑线表示。图中红线是卡尔曼滤波器输出的最优化结果（该结果在算法中设置了Q=1e-6，R=1e-1）。第三节卡尔曼滤波的信号模型（Signal Model of

50、 Kalman Filtering） 通过前面几节内容的学习，我们知道维纳滤波是根据当前 和过去全部的观测值 来估计信号的当前值 ，它的解形式是以均方误差最小为原则下的系统的传递函数 或单位脉冲响应 。 卡尔曼滤波不需要过去全部的观测, 它是根据前一个估计值)(nx),2(),1(nxnx)(ns)(zH)(nh) 1( ns)(nx)( ns和最近一个观测值来估计信号的当前推方法进行估计的，因而卡尔曼滤波对信号的平稳性和时不变性不做要求。我们利用维纳滤波的模型引入到卡尔曼滤波的信号模型。，它是用状态方程和递一、状态一、状态方程和量方程和量测方程测方程 要给出卡尔曼滤波的信号模型，先来讨论状态

51、方程和量测方程。图9-11是维纳滤波的模型，信号 可以认为是由白噪声 激励一个线性系统 的响应，假设响应和激励的时域关系可以用下式表示： (9-52) 1() 1()(1nwnasns)(ns)(1nw)(zA上式也就是一阶AR模型。在卡尔曼滤波中信号)(ns被称为是状态变量，用矢量的 S(k)S(k)形式表示为，在k时刻的状态用表示， 在k1时刻的状态用1)S(k 表示。 )(ns)(zA)(1nw)(ns)(zA)(1nw)(nw)(nx图9-11 维纳滤波的信号模型和观测信号模型)(1nw(k)w1A(k)激励信号也用矢量表示为，激励和响应之间的关系用传递矩阵来表 示，它是由系统的)(z

52、A有一定关系。有了这些假设后结构确定的，与我们给出状态方程：1)(kw1)A(k)S(kS(k)1 (9-53)S(k)1)S(k 1)S(k 上式表示的含义就是在k时刻的状态可以由它的前一个时刻的状态来求得，即认为k1时刻以前的各状态都已记忆在状态中了 )(ns)(zA)(1nw)(ns)(zA)(1nw)(nw)(nx图9-11 维纳滤波的信号模型和观测信号模型)()()(nwnsnx)(nw卡尔曼滤波是根据系统的量测数据（即观测数据）对系统的运动进行估计的，所以除了状态方程之外，还需要量测方程。还是从维纳滤波的观测信号模型入手，图9-11的右图，观测数据和信号的关系为：，一般是均值为零的

53、高斯白噪声，卡尔曼滤波中X(k)S(k)，量测矢量与状态矢量的关系为w(k)S(k)X(k) (9-54)上式和维纳滤波的)()()(nwnsnx概念上是一致的，也就是说卡尔曼滤波的一维信号模型和维纳滤波的信号模型是一致的。把式(9-54)推广就得到更普遍的多维量测方程w(k)C(k)S(k)X(k) (9-55)C(k)X(k)S(k)上式中的称为量测矩阵，它的引入原因是，的维数不一定与状态矢量的维数相同，因为我们不一定能观测到所有需要的状态参数 量测矢量X(k)1mS(k)1nC(k)nmw(k)1m假如是的矢量，是的矢量，就是的矩阵，是的矢量。二、信号模型 有了状态方程 和量测方程 后我

54、们就能给出卡尔曼滤波的信号模型，如图9-12所示。1)(kw1)A(k)S(kS(k)1w(k)C(k)S(k)X(k)S(k)C(k)1)A(k 1zw(k)(k)w1X(k)1)S(k 图9-12 卡尔曼滤波的信号模型w(k)S(k)X(k)25. 18 . 0 ,)8 . 01)(8 . 01 (36. 0)(1zzzzRss)()(mmRwwA(k)C(k)【例9-6】设卡尔曼滤波中量测方程为，已知信号的自相关函数的z变换为 噪声的自相关函数为：，信号和噪声统计独立。求卡尔曼滤波和。信号模型中的)()()(121zAzAzRwss)8 . 01)(8 . 01 (36. 011zzzz

55、解根据等式：可以求得：)()(8 . 01)(111zWzSzzzA)()(8 . 0) 1(1kwksks8 . 0A(k)w(k)S(k)X(k)C(k)变换到时域得：因此 又因为，所以1。第五节第五节 卡尔曼滤波方法卡尔曼滤波方法（Method of Kalman Filtering） 建立好了卡尔曼滤波的信号模型以及状态方程、量测方程后，要解决的问题就是要寻找在最小均方误差下信号 的估计值 。S(k)(k)S1)(kw1)A(k)S(kS(k)1 (9-56)w(k)C(k)S(k)X(k) (9-57)A(k)C(k)X(k)1)(kS(k)S上式中和是已知的，已知，现在的问题就是如

56、何来求当前时刻。 是观测到的数据，也是已知的，假设信号的 上一个估计值 的估计值(k)w1w(k)S(k)(k)w1w(k)(k)S(k)X(k)S(k)X上两式中如果没有与，可以立即求得，估计问题的出现就是因为信号与噪声的与，用上两式和分别用和表示，得：叠加。假设暂不考虑得到的1)(kSA(k)(k)S (9-58)1)(kSC(k)A(k)(k)SC(k)(k)X (9-59)X(k)(k)X必然，观测值和估计值之间有误差 ，它们之间的差(k)X称为新息（innovation）：(k)XX(k)(k)X (9-60)(k)w1w(k)显然，新息的产生是由于我们前面忽略了与所引起的，也就是说

57、新息里面 (k)w1w(k)(k)XH(k)(k)w1S(k)包含了与的信息成分。因而我们用新息乘以一个修正矩阵，用它来代替式（9-56）的来对进行估计：(k)XH(k)1)(kSA(k)(k)S1)(kSC(k)A(k)H(K)X(k)1)(kSA(k) (9-61)S(k)X(k)(k)S(k)SC(k)A(k)1zX(k)1)(kS H(k)(k)X (k)X由（9-56）（9-61）可以画出卡尔曼滤波对进行估计的递推模型，如图9-13所示，输入为观测值，输出为信号估计值。图9-13 卡尔曼滤波的一步递推法模型二、卡尔曼滤波的递推公式 从图9-13容易看出，要估计出 就必须要先找到最小均

58、方误差下的修正矩阵 ，结合式（9-61）、（9-56）、（9-57）得：(k)SH(k)1)(kSC(k)A(k)w(k)(k)H(K)C(k)S1)(kSA(k)(k)S1)(kSC(k)A(k)w(k)1)(kw1)(kSA(k)H(K)C(k)1)(kSA(k)1H(k)w(k)1)(kw1)(kS(k)H(K)C(k)AH(k)C(k)1)I(kSA(k)1.(9-62)H(k)H(k)根据上式来求最小均方误差下的，然后把求到的代入（9-61）则可以得 到估计值(k)S。(k)SS(k)(k)S(k)设真值和估计值之间的误差为：，误差是个矢量，因而均方误差是一个矩阵，用表示。把式（9-

59、62）代入得：(k)SS(k)(k)SH(k)w(k)1)(kw1)(kS1)A(k)S(kH(K)C(k)I1.(9-63)均方误差矩阵：(k)S(k)SE(k) (9-64)表示对向量取共轭转置。为了计算方便，令(k)S(k)(S(k)SE(S(k)(k) (9-65) 找到和均方误差矩阵的关系：1)(kSA(k)1)(kw1)k1)(A(k)S(kSA(k)1)(kw1)E(A(k)S(k(k)111)(k1)w(kEwA(k)1)(kS1)1)(S(k(kS1)A(k)ES(k111)Q(k1)A(k)A(k)(k (9-66)(k)w1w(k)j)Q(k)(k(j)(k)wEw11把

60、式（9-63）代入式（9-64），并且利用条件：与都是零均值的高斯白噪声，且它们和 互不相关，协方差矩阵分别为j)R(k)(kEw(k)w(j)1)S(k)(kw111)(kS)(kw11w(k)与不相关；与及不相关。最后化简得：(k)S(k)SE(k)k)H(k)R(k)H(H(k)C(k)1)IQ(k1)A(k)A(k)(kH(K)C(k)I. (9-67)把式（9-66）代入（9-67）得(k)k)H(k)R(k)H(H(k)C(k)(k)IH(K)C(k)IR(k)H(k)(k)C(k)H(k)C(k)H(k)(k)C(k)(k)H(K)C(k)(k)SSR(k)(k)C(k)C(k)