《东南大学统计信号处理实验实验三.docx》由会员分享，可在线阅读，更多相关《东南大学统计信号处理实验实验三.docx（19页珍藏版）》请在三一文库上搜索。
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); X0
5、(:,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,X(2,:),b,t,X0(2,:),r);h=legend(真实值,估计值);xlabel(t);ylabel(v1);title(alfa-be
6、ta滤波v1);grid on;figure(3);plot(t,X(4,:),b,t,X0(4,:),r);h=legend(真实值,估计值);xlabel(t);ylabel(v2);title(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,4
7、99);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 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
8、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=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
9、)=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)=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,:)
10、,X0(1,:),r);h=legend(真实值,估计值);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,X0(4,:),r);h=legend(真实值,估计值);xlabel(t);ylabel(v2);title(Kalman滤波的v2);grid on;3)假设目标在运动过程
11、中发生了机动（速度在某个时刻突然发生了改变），试观测此时的滤波和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 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=
12、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*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)+10; X(4,i+1)= X(4,i+1)+10; end Y(
13、:,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); 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=le
14、gend(真实值,估计值);xlabel(s2);ylabel(s1);title(alfa-beta滤波的距离);grid on;figure(2);plot(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);title(alfa-beta滤波的v2);grid on;卡尔曼滤波clear;clc;sig
15、ma_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 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
16、;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=sigma_u1*randn(1,500);u2=sigma_u2*randn(1,500);w1=sigma_w*randn(1
17、,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(
18、:,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(真实值,估计值);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,X0(4,:),r);h=legend(真实值,估计值);xlabel(t);ylabel(v2);title(Kalman滤波的v2);grid on;