东南大学统计信号处理实验实验三

1、 统计信号处理实验三目的：掌握卡尔曼滤波滤波器的原理；内容：用雷达跟踪目标，目标的运动可以看成是在径向和横向内的二维运动，其运动方程和观测方程分别为：、和分别为径向距离、速度和观测值，而、和分别为横向距离、速度和观测值。和是状态噪声，是目标速度的波动；和是观测噪声；四种噪声的均值都为0，呈高斯分布，互不相关。T是雷达扫描一次的时间，此处设为1.0秒。假设目标距离雷达约160Km左右，径向初速度设为300 m/s，并且在向雷达靠近，横向初速度设为0 m/s。这样它的径向速度波动大，而横向速度波动小，所以我们假设的方差为300m/s，的方差为m/s。鉴于雷达的观测误差，我们假设观测噪声和的方差和均

2、为1.0Km。其中，和的初始值不是最佳的，学生完全可自己修改以上参数，并观察计算结果的变化，给出最好的滤波效果。任务：1） 试用滤波法对信号进行处理，并通过计算机模拟对其跟踪过程进行验证；2） 试求其Kalman滤波方程，并通过计算机模拟对其跟踪过程进行验证；3） 假设目标在运动过程中发生了机动（速度在某个时刻突然发生了改变），试观测此时的滤波和Kalman滤波结果，并对结果进行解释。要求：1）设计仿真计算的Matlab程序，给出软件清单；2）完成实验报告，给出实验结果，并对实验数据进行分析。(1) 滤波法对信号进行处理。clear;alfa=0.6;beta=0.4;sigma_u1=0.3

3、;sigma_u2=0.2;sigma_w=0.1;T=1; s1=12 zeros(1,499);s2=7 zeros(1,499);v1=15 zeros(1,499);v2=4 zeros(1,499);y1=zeros(1,500);y2=zeros(1,500);A=1 T 0 0;0 1 0 0;0 0 1 T;0 0 0 1;C=1 0 0 0;0 0 1 0;X=s1;v1;s2;v2;X0=11.8 zeros(1,499);13.8 zeros(1,499);6.8 zeros(1,499);3.9 zeros(1,499);Y=y1;y2;Y0=y1;y2;u1=sigm

4、a_u1*randn(1,500);u2=sigma_u2*randn(1,500);w1=sigma_w*randn(1,500);w2=sigma_w*randn(1,500);M=500;for i=1:M-1; X(:,i+1)=A*X(:,i)+0;u1(i);0;u2(i); Y(:,i+1)=C*X(:,i+1)+w1(i+1) w2(i+1)'ends1=X(1,:);v1=X(2,:);s2=X(3,:);v2=X(4,:);K=alfa 0;beta/T 0;0 alfa;0 beta/T;M=500;for i=1:M-1; X1(:,i+1)=A*X0(:,i)

5、; X0(:,i+1)=X1(:,i+1)+K*Y(1,i+1)-X1(1,i+1);Y(2,i+1)-X1(3,i+1); Y0(:,i+1)=C*X0(:,i+1);endt=0:499;figure(1);plot(X(3,:),X(1,:),'b',X0(3,:),X0(1,:),'r');h=legend('真实值','估计值');xlabel('s2');ylabel('s1');title('alfa-beta滤波');grid on;figure(2);plot(t,

6、X(2,:),'b',t,X0(2,:),'r');h=legend('真实值','估计值');xlabel('t');ylabel('v1');title('alfa-beta滤波v1');grid on;figure(3);plot(t,X(4,:),'b',t,X0(4,:),'r');h=legend('真实值','估计值');xlabel('t');ylabel('v2');ti

7、tle('alfa-beta滤波v2');grid on;取值alfa = 0.8；beta = 0.2；第二个取值下的估计效果较差。2)试求其Kalman滤波方程，并通过计算机模拟对其跟踪过程进行验证；clear;clc;sigma_u1=0.3;sigma_u2=0.2;sigma_w=0.1;T=1;s1=12 zeros(1,499);s2=7 zeros(1,499);v1=15 zeros(1,499);v2=4 zeros(1,499);y1=zeros(1,500);y2=zeros(1,500);A=1 T 0 0;0 1 0 0;0 0 1 T;0 0 0

8、1;C=1 0 0 0;0 0 1 0;Q=0 0 0 0;0 sigma_u1 0 0;0 0 0 0;0 0 0 sigma_u2;R=sigma_w 0;0 sigma_w;I=1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1;P2=sigma_w sigma_w/T 0 0; sigma_w/T sigma_u1+2*sigma_w/(T2) 0 0; 0 0 sigma_w sigma_w/T; 0 0 sigma_w/T sigma_u2+2*sigma_w/(T2) ;X=s1;v1;s2;v2;X0=11.8 zeros(1,499);13.8 zeros(1,4

9、99);6.8 zeros(1,499);3.9 zeros(1,499);Y=y1;y2;Y0=y1;y2;u1=sigma_u1*randn(1,500);u2=sigma_u2*randn(1,500);w1=sigma_w*randn(1,500);w2=sigma_w*randn(1,500);M=500;for i=1:M-1; X(:,i+1)=A*X(:,i)+0;u1(i);0;u2(i); Y(:,i+1)=C*X(:,i+1)+w1(i+1) w2(i+1)'ends1=X(1,:);v1=X(2,:);s2=X(3,:);v2=X(4,:);for i=1:M-

10、1; if i=1 P1=A*P2*A'+Q; else P1=A*P0*A'+Q; end K=P1*C'*inv(C*P1*C'+R); X1(:,i+1)=A*X0(:,i) X0(:,i+1)=X1(:,i+1)+K*(Y(:,i+1)-C*X1(:,i+1); Y0(:,i+1)=C*X0(:,i+1); P0=(I-K*C)*P1;endt=1:500;figure(1);plot(X(3,:),X(1,:),'b',X0(3,:),X0(1,:),'r');h=legend('真实值','估计

11、值');xlabel('s2');ylabel('s1');title('Kalman滤波的距离');grid on;figure(2);plot(t,X(2,:),'b',t,X0(2,:),'r');h=legend('真实值','估计值');xlabel('t');ylabel('v1');title('Kalman滤波的v1');grid on;figure(3);plot(t,X(4,:),'b',t,

12、X0(4,:),'r');h=legend('真实值','估计值');xlabel('t');ylabel('v2');title('Kalman滤波的v2');grid on;3)假设目标在运动过程中发生了机动（速度在某个时刻突然发生了改变），试观测此时的滤波和Kalman滤波结果，并对结果进行解释。滤波clear;clc;alfa=0.6;beta=0.4;sigma_u1=0.3;sigma_u2=0.2;sigma_w=0.1;T=1; s1=12 zeros(1,499);s2=7 zer

13、os(1,499);v1=15 zeros(1,499);v2=4 zeros(1,499);y1=zeros(1,500);y2=zeros(1,500);A=1 T 0 0;0 1 0 0;0 0 1 T;0 0 0 1;C=1 0 0 0;0 0 1 0;X=s1;v1;s2;v2;X0=11.8 zeros(1,499);13.8 zeros(1,499);6.8 zeros(1,499);3.9 zeros(1,499);Y=y1;y2;Y0=y1;y2;u1=sigma_u1*randn(1,500);u2=sigma_u2*randn(1,500);w1=sigma_w*rand

14、n(1,500);w2=sigma_w*randn(1,500);M=500;for i=1:M-1; X(:,i+1)=A*X(:,i)+0;u1(i);0;u2(i); if i=199 X(2,i+1)= X(2,i+1)+10; X(4,i+1)= X(4,i+1)+10; end Y(:,i+1)=C*X(:,i+1)+w1(i+1) w2(i+1)' ends1=X(1,:);v1=X(2,:);s2=X(3,:);v2=X(4,:);K=alfa 0;beta/T 0;0 alfa;0 beta/T;M=500;for i=1:M-1; X1(:,i+1)=A*X0(:,

15、i); X0(:,i+1)=X1(:,i+1)+K*Y(1,i+1)-X1(1,i+1);Y(2,i+1)-X1(3,i+1); Y0(:,i+1)=C*X0(:,i+1);endt=0:499;figure(1);plot(X(3,:),X(1,:),'b',X0(3,:),X0(1,:),'r');h=legend('真实值','估计值');xlabel('s2');ylabel('s1');title('alfa-beta滤波的距离');grid on;figure(2);pl

16、ot(t,X(2,:),'b',t,X0(2,:),'r');h=legend('真实值','估计值');xlabel('t');ylabel('v1');title('alfa-beta滤波的v1');grid on;figure(3);plot(t,X(4,:),'b',t,X0(4,:),'r');h=legend('真实值','估计值');xlabel('t');ylabel('v2

17、9;);title('alfa-beta滤波的v2');grid on;卡尔曼滤波clear;clc;sigma_u1=0.3;sigma_u2=0.2;sigma_w=0.1;T=1;s1=12 zeros(1,499);s2=7 zeros(1,499);v1=15 zeros(1,499);v2=4 zeros(1,499);y1=zeros(1,500);y2=zeros(1,500);A=1 T 0 0;0 1 0 0;0 0 1 T;0 0 0 1;C=1 0 0 0;0 0 1 0;Q=0 0 0 0;0 sigma_u1 0 0;0 0 0 0;0 0 0 si

18、gma_u2;R=sigma_w 0;0 sigma_w;I=1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1;P2=sigma_w sigma_w/T 0 0; sigma_w/T sigma_u1+2*sigma_w/(T2) 0 0; 0 0 sigma_w sigma_w/T; 0 0 sigma_w/T sigma_u2+2*sigma_w/(T2) ;X=s1;v1;s2;v2;X0=11.8 zeros(1,499);13.8 zeros(1,499);6.8 zeros(1,499);3.9 zeros(1,499);Y=y1;y2;Y0=y1;y2;u1=sig

19、ma_u1*randn(1,500);u2=sigma_u2*randn(1,500);w1=sigma_w*randn(1,500);w2=sigma_w*randn(1,500);M=500;for i=1:M-1; X(:,i+1)=A*X(:,i)+0;u1(i);0;u2(i); if i=199 X(2,i+1)= X(2,i+1)+8; X(4,i+1)= X(4,i+1)+8; end Y(:,i+1)=C*X(:,i+1)+w1(i+1) w2(i+1)'ends1=X(1,:);v1=X(2,:);s2=X(3,:);v2=X(4,:);for i=1:M-1; if i=1 P1=A*P2*A'+Q; else P1=A*P0*A'+Q; end K=P1*C'*inv(C*P1*C'+R); X1(:,i+1)