数字信号处理实验2

数字信号处理实验 实验二、Kalman滤波估计

一、实验内容及要求 如果kkkkkkvxywxx1, 试求',,,ˆkkkkPPHx(k=1,2,…,6,7)，并画出曲线。

其中w和v是零均值不相关过程，且有kjjkwwCov],[，kjjkvvCov],[，00，观测序列}1,1,1,3,4,2,1{)7(y，并考虑以下三种情况 (1) 0P, (2) 10P , (3) 00P 二、问题解决思路 波的信号模型是从状态方程和量测方程得到。设动态系统的状态方程和量测方程分别为

11kkkkxAx (1)

kkkkxCy (2)

其中： kA——nn维矩阵，状态转移矩阵；

kC——nm维矩阵，量测矩阵；

kx——n维状态向量；

ky——m维观测向量；

kw——n维（过程）噪声向量；

kv——m维（量测）噪声向量；

假设kk与都是均值为零的正态白噪声，且kk与互不相关，即 kjkjkjkkQEE , ,cov ;0 (3)

kjkjkjkkREE , ,cov ;0 (4)

 2, ,1 ,0 , ;0 , ,cov jkEjkjk

这里 

kkkkkkkkEREQvar ;var

并设初始状态0x与jk 、均不相关。 上面这些假设是符合一般实际情况的，根据这样的假设，在满足均方误差阵最小的条件下，我们得到下列一组卡尔曼一步递推公式：

)ˆ(ˆˆ11kkkkkkkkxACyHxAx (5)

1)(kkkkkkkRCPCCPH



(6)

11kkkkkQAPAP (7)

kkkkPCHIP (8)

其中： kH为利用前一个估计值和当前观测数据对预测估计值进行修正的修正矩阵；



kkkkkkkxxExxxxEP~ ~)ˆ)(ˆ( (9)

并令 )ˆ)(ˆ(kkkkkxxxxEP

(10)

1ˆˆkkkxAx (11)

由式(5)可见，如果已知Hk ，利用前一个kx的估计值1ˆkx与当前的量测值ky，就可以求得kxˆ。如果Hk是按式(6)计算的，即满足最小均方误差阵的Hk ，则将此Hk代入式(5)，就得

到我们所要求的在最小均方误差阵条件下的kxˆ。如果初始状态x0的系统特性00varxxE及已知，并

令 000ˆxEx

又 000000varˆˆxxxxxEP 将P0代入式(7)可求得1P，将1P代入式(6)可求得H1，将此H1代入式(5)可求得在最小均方误差条件下的1ˆx，同时将1P代入式(8)又可求得P1；由P1又可求得2P，由2P又可求得H2，由H2又可求得2ˆx，同时由H2与2P又可求得P2……；以此类推。故根据已知的矩阵：kkkkkyRQCA以及观测到的、、、，就能用递推计算法得到所有的kxxxˆ ˆˆ21，、、以及

kPPP、、、21。

根据题目中给出的状态方程和量测方程可以得到kkkkRQCA、、、的值分别为： 1kA，1kC，1kQ，1kR 由此有给定信号模型下的Kalman滤波一步递推公式：







kkkkkkkkkkkkkPHPPPPPHxyHxx)1(

1)1()ˆ(ˆˆ1111

通过代入消去kP，得到kkPH，同时可以求得Kalman过滤达到稳态时 618.025112114112P 即，有kP将逐步趋向于0.618。将初值0ˆ00x，观测序列}1,1,1,3,4,2,1{)7(y及(1)0P, (2) 10P, (3)00P分别代入Kalman一步递推公式即可得到所求的',,,ˆkkkkPPHx(k=1,2,…,6,7)。

三、实验步骤及结果分析 a) 流程图

(12) (13) (14) (15) b) 令0P分别取∞、1和0，}1,1,1,3,4,2,1{)7(y，计算',,,ˆkkkkPPHx(k=1,2,…,6,7)，并画出它们的变化曲线，如图1~4所示。所求得的',,,ˆkkkkPPHx(k=1,2,…,6,7)为： （1）0P=∞

kxˆ={1.0000 -1.0000 2.1250 2.6667 0.4000 0.7708 0.9125}

kH={1.0000 0.6667 0.6250 0.6190 0.6182 0.6181 0.6180}

kP={1.0000 0.6667 0.6250 0.6190 0.6182 0.6181 0.6180}

'kP={ ∞ 2.0000 1.6667 1.6250 1.6190 1.6182 1.6181}

（2）0P=1 kxˆ={0.6667 -1.0000 2.0952 2.6545 0.3958 0.7692 0.9119}

kH={0.6667 0.6250 0.6190 0.6182 0.6181 0.6180 0.6180}

kP={0.6667 0.6250 0.6190 0.6182 0.6181 0.6180 0.6180}

变量定义赋初值 P0取∞，1，0 K = 1~7 P0 ? ∞ 将数据代入卡尔曼一步递推公式

结束 描点画图

H1=1 取遍 假 未取遍

真 是 否 'kP={2.0000 1.6667 1.6250 1.6190 1.6182 1.6181 1.6180}

（3）0P=0 kxˆ={0.5000 -1.0000 2.0769 2.6471 0.3933 0.7682 0.9115}

kH={0.5000 0.6000 0.6154 0.6176 0.6180 0.6180 0.6180}

kP={0.5000 0.6000 0.6154 0.6176 0.6180 0.6180 0.6180}

'kP={1.0000 1.5000 1.6000 1.6154 1.6176 1.6180 1.6180}

实验结果分析：

(1) 实验结果表明，在0P分别取∞、1和0的情况下，Kalman估计结果随迭代次数的增多而趋于同一值，kP逐步趋于P。可见，Kalman过滤有一个过渡过程，对于线性时不变系统，Kalman滤波只有达到稳态时才与维纳滤波有相同的结果，且稳态下的方均误差P与其初值00varxP无关。

(2) 从所画出的曲线可以观察到，kH、kP和'kP都存在一个极限值，理论计算公式推导得到'kP≥Q(=1)，0≤kH≤1，0≤kP≤R(=1)，即，jk 、方差的大小分别决定了预测精度的下限和滤波估计精度的上限。显然，实验结果都是满足的。k

H

的变化趋势表明，当预测估计误差逐渐减小时，修正参数也逐渐减小，直到达到一个估计系统无法逾越的极限。 图1、0P分别取∞、1和0时，kxˆ（k=1,2…6,7）的变化曲线

图2、0P分别取∞、1和0时，修正参数kH（k=1,2…6,7）的变化曲线 图3、0P分别取∞、1和0时，kP（k=1,2…6,7）的变化曲线

图4、0P分别取∞、1和0时，'kP（k=1,2…6,7）的变化曲线 四、 源程序 %% dsp homework2 Kalman filter design %% designed by HanBing 200328013703525 clear; close all; Ex = zeros(8,1); % estimation of expect signal y = [1;-2;4;3;-1;1;1]; % observed signal H = zeros(7,1); % adjusting parameter P = zeros(8,1); % square error of Ex and x Pe = zeros(7,1); % square error of Ex' and x temp = [inf , 1 , 0]; colr = ['r','+',':' ; 'g','+',':' ; 'b','+',':']; for i = 1:3 P(1) = temp(i); for k = 2:8 Pe(k-1) = P(k-1) + 1; if Pe(k-1) == inf H(k-1) = 1; P(k) = 1; else H(k-1) = Pe(k-1) / (Pe(k-1)+1); P(k) = H(k-1); end Ex(k) = Ex(k-1) + H(k-1) * (y(k-1) - Ex(k-1)); end t = 0:7; figure(1); plot(t,Ex,colr(i,:)); hold on; title('Ex(k) ( P0 = [ Inf (red) ; 1 (green) ; 0 (blue) ] )');

figure(2); plot(H,colr(i,:)); hold on; title('H(k) ( P0 = [ Inf (red) ; 1 (green) ; 0 (blue) ] )');

figure(3); plot(t,P,colr(i,:)); hold on; title('P(k) ( P0 = [ Inf (red) ; 1 (green) ; 0 (blue) ] )');

figure(4); plot(Pe,colr(i,:)); hold on;