东南大学统计信号处理实验实验三
东南大学统计信号处理实验一
《统计信号处理》实验一一、实验目的：1、掌握噪声中信号检测的方法；2、熟悉Matlab 的使用；3、掌握用计算机进行数据分析的方法。
二、实验内容：假设信号为()s t 波形如下图所示：在有信号到达时接收到的信号为()()()x t s t n t =+，在没有信号到达时接收到的信号为()()x t n t =。
其中()n t 是均值为零、方差为225n σ=（可自行调整）的高斯白噪声。
假设有信号到达的概率P(H 1)=0.6，没有信号到达的概率P(H 0)=0.4。
对接受到的信号分别在t = 0ms, 1ms, …, 301ms 上进行取样，得到观测序列()x n 。
1、利用似然比检测方法（最小错误概率准则），对信号是否到达进行检测；2、假设102C =，011C =。
利用基于Bayes 准则的检测方法，对信号是否到达进行检测；3、通过计算机产生的仿真数据，对两种方法的检测概率d P 、虚警概率f P 、漏警概率m P 和Bayes 风险进行仿真计算；4、通过改变P(H 1)和P(H 0)来改变判决的门限（风险系数10C 和01C 不变），观察检测方法的d P 、f P 、m P 和Bayes 风险的变化；5、改变噪声的方差，观察检测方法的d P 、f P 、m P 和Bayes 风险的变化；6、将信号取样间隔减小一倍(相应的取样点数增加一倍)，观察似然比检测方法的d P 、f P 、m P 和Bayes 风险的变化；7、根据()s t 设计一个离散匹配滤波器，并观察()x n 经过该滤波器以后的输出。
三、实验要求：1、设计仿真计算的Matlab 程序，给出软件清单；2、完成实验报告，对实验过程进行描述，并给出实验结果，对实验数据进行分析，给出结论。
四、设计过程：1、产生信号s(t)，n(t)，x(t)，t = 0ms, 1ms, …, 301ms ；其中：⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧≤≤+-≤≤-≤≤+-≤≤-≤≤+-≤≤-≤≤+-≤≤=301290,30101289270,28101269230,5.12201229190,5.10201189140,6.625113990,6.42518930,2301290,301)(t t t t t t t t t t t t t t t t t s 2、根据定义似然比函数10(|)()(|)p x H x p x H Λ=，门限001()()P H P H Λ=，如果0)(Λ>Λx ，则判定1D ；否则，判定0D 。
东南大学检测实验报告
传感器第一次实验试验一 金属箔式应变片——单臂电桥性能实验一． 实验目的了解金属箔式应变片的应变效应及单臂电桥工作原理和性能。
二． 基本原理电阻丝在外力作用下发生机械形变时，其电阻值发生变化，这就是电阻应变效应。
金属箔式应变片就是通过光刻、腐蚀等工艺制成的应变敏感元件，通过它反映被测部位受力状态的变化。
电桥的作用是完成电阻到电压的比例变化，电桥的输出电压反映了相应的受力状态。
单臂电桥输出电压 1/4o U EK ε=，其中K 为应变灵敏系数，/L L ε=∆为电阻丝长度相对变化。
三． 实验器材主机箱、应变传感器实验模板、托盘、砝码、万用表、导线等。
四． 实验步骤1. 根据接线示意图安装接线。
2. 放大器输出调零。
3. 电桥调零。
4. 应变片单臂电桥实验。
测得数据如下： 重量（g ） 0 20 40 60 80 100 120 140 电压（mv ）3.57.211.316.020.724.728.833.5实验曲线如下所示：分析：由图可以看出，输出电压与加载的重量成线性关系，由于一开始调零不好，致使曲线没有经过原点，往上偏离了一段距离。
5. 根据表中数据计算系统的灵敏度/S U W =∆∆（U ∆为输出电压变化量，W ∆为重量变化量）和非线性误差/100%m yFS δ=∆⨯，式中m ∆为输出值（多次测量时为平均值）与拟合直线的最大偏差；yFS 为满量程输出平均值，此处为140g 。
U ∆=30mv ， W ∆=140g ， 所以 30/1400.2143/S mv g == m ∆=1.9768g ， yFS =140g ， 所以 1.9768/140100% 1.41%δ=⨯=6. 利用虚拟仪器进行测量。
测得数据如下表所示： 重量（g ） 0 20 40 60 80 100 120 140 电压（mv ）0.75.09.513.918.723.428.332.9相应的曲线如下：五． 思考题单臂电桥工作时，作为桥臂电阻的应变片应选用：（1）正（受拉）应变片；（2）负（受压）应变片；（3）正、负应变片均可以。
数字信号处理实验
数字信号处理实验实验一信号、系统及系统响应1、实验目的认真复习采样理论、离散信号与系统、线性卷积、序列的z 变换及性质等有关内容；掌握离散时间序列的产生与基本运算，理解离散时间系统的时域特性与差分方程的求解方法，掌握离散信号的绘图方法；熟悉序列的z 变换及性质，理解理想采样前后信号频谱的变化。
2、实验内容a. 产生长度为500 的在[0,1]之间均匀分布的随机序列，产生长度为500 的均值为0 单位方差的高斯分布序列。
b. 线性时不变系统单位脉冲响应为h(n)=(0.9)nu(n)，当系统输入为x(n)=R10(n)时，求系统的零状态响应，并绘制波形图。
c. 描述系统的差分方程为：y(n)-y(n-1)+0.9y(n-2)=x(n)，其中x(n)为激励，y(n)为响应。
计算并绘制n=20,30,40,50,60,70,80,90,100 时的系统单位脉冲响应h(n)；计算并绘制n=20,30,40,50,60,70,80,90,100 时的系统单位阶跃响应s(n)；由h(n)表征的这个系统是稳定系统吗？d. 序列x(n)=(0.8)nu(n)，求DTFT[x(n)]，并画出它幅度、相位，实部、虚部的波形图。
观察它是否具有周期性？e. 线性时不变系统的差分方程为y(n)=0.7y(n-1)+x(n)，求系统的频率响应H(ejω)，如果系统输入为x(n)=cos(0.05πn)u(n)，求系统的稳态响应并绘图。
f. 设连续时间信号x(t)=e-1000|t|，计算并绘制它的傅立叶变换；如果用采样频率为每秒5000 样本对x(t)进行采样得到x1(n)，计算并绘制X1(ejω)，用x1(n)重建连续信号x(t)，并对结果进行讨论；如果用采样频率为每秒1000 样本对x(t)进行采样得到x2(n)，计算并绘制X2(ejω)，用x2(n)重建连续信号x(t)，并对结果进行讨论。
加深对采样定理的理解。
g. 设X1(z)=z+2+3z-1，X2(z)=2z2+4z+3+5z-1，用卷积方法计算X1(z)X2(z)。
东南大学信号与系统MATLAB实践第三次作业
练习三实验三五．1．>>help windowWINDOW Window function gateway.WINDOW(@WNAME,N) returns an N-point window of type specifiedby the function handle @WNAME in a column vector. @WNAME canbe any valid window function name, for example:@bartlett - Bartlett window.@barthannwin - Modified Bartlett-Hanning window.@blackman - Blackman window.@blackmanharris - Minimum 4-term Blackman-Harris window.@bohmanwin - Bohman window.@chebwin - Chebyshev window.@flattopwin - Flat Top window.@gausswin - Gaussian window.@hamming - Hamming window.@hann - Hann window.@kaiser - Kaiser window.@nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window.@parzenwin - Parzen (de la Valle-Poussin) window.@rectwin - Rectangular window.@tukeywin - Tukey window.@triang - Triangular window.WINDOW(@WNAME,N,OPT) designs the window with the optional input argument specified in OPT. To see what the optional input arguments are, see the helpfor the individual windows, for example, KAISER or CHEBWIN.WINDOW launches the Window Design & Analysis Tool (WinTool).EXAMPLE:N = 65;w = window(@blackmanharris,N);w1 = window(@hamming,N);w2 = window(@gausswin,N,2.5);plot(1:N,[w,w1,w2]); axis([1 N 0 1]);legend('Blackman-Harris','Hamming','Gaussian');See also bartlett, barthannwin, blackman, blackmanharris, bohmanwin,chebwin, gausswin, hamming, hann, kaiser, nuttallwin, parzenwin,rectwin, triang, tukeywin, wintool.Overloaded functions or methods (ones with the same name in other directories) help fdesign/window.mReference page in Help browser doc window 2.>>N = 128;w = window(@rectwin,N); w1 = window(@bartlett,N); w2 = window(@hamming,N);plot(1:N,[w,w1,w2]); axis([1 N 0 1]); legend('矩形窗','Bartlett','Hamming');2040608010012000.10.20.30.40.50.60.70.80.913.>>wvtool(w,w1,w2)六．ts=0.01;N=20;t=0:ts:(N-1)*ts;x=2*sin(4*pi*t)+5*cos(6*pi*t);g=fft(x,N);y=abs(g)/100;figure(1):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid;01234560.050.10.150.20.250.30.350.40.45ts=0.01; N=30;t=0:ts:(N-1)*ts;x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100;figure(2):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid;012345670.10.20.30.40.50.60.7ts=0.01; N=50;t=0:ts:(N-1)*ts;x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100;figure(3):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid;123456700.10.20.30.40.50.60.70.80.91ts=0.01; N=100;t=0:ts:(N-1)*ts;x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100;figure(4):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid;012345670.511.522.5ts=0.01; N=150;t=0:ts:(N-1)*ts;x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100;figure(5):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid;012345670.511.522.53实验八 1．%冲激响应 >> clear; b=[1,3]; a=[1,3,2]; sys=tf(b,a); impulse(sys); 结果：Impulse ResponseTime (sec)A m p l i t u d e%求零输入响应 >> A=[1,3;0,-2]; B=[1;2]; Q=A\B Q = 4-1>> clear B=[1,3]; A=[1,3,2];[a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[4;-1]; initial(sys,x0); grid; a =-3 -2 1 0 b =1 0 c =1 3 d = 001234560.20.40.60.811.21.4Response to Initial ConditionsTime (sec)A m p l i t u d e2.%冲激响应 >> clear; b=[1,3]; a=[1,2,2]; sys=tf(b,a); impulse(sys)Impulse ResponseTime (sec)A m p l i t u d e%求零输入响应 >> A=[1,3;1,-2]; B=[1;2]; Q=A\B Q =1.6000 -0.2000 >> clear B=[1,3]; A=[1,2,2];[a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[1.6;-0.2]; initial(sys,x0); grid; a =-2 -2 1 0 b = 10 c =1 3 d = 0Response to Initial ConditionsTime (sec)A m p l i t u d e3.%冲激响应 >> clear; b=[1,3]; a=[1,2,1]; sys=tf(b,a); impulse(sys)Impulse ResponseTime (sec)A m p l i t u d e%求零输入响应 >> A=[1,3;1,-1]; B=[1;2]; Q=A\B Q =1.7500 -0.2500>> clear B=[1,3]; A=[1,2,1];[a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[1.75;-0.25]; initial(sys,x0); grid; a =-2 -1 1 0 b =1 0 c =1 3 d = 00510150.20.40.60.811.21.41.6Response to Initial ConditionsTime (sec)A m p l i t u d e二． >> clear; b=1;a=[1,1,1,0]; sys=tf(b,a); subplot(2,1,1);impulse(sys);title('冲击响应'); subplot(2,1,2);step(sys);title('阶跃响应'); t=0:0.01:20; e=sin(t);r=lsim(sys,e,t); figure;subplot(2,1,1);plot(t,e);xlabel('Time');ylabel('A');title('激励信号'); subplot(2,1,2);plot(t,r);xlabel('Time');ylabel('A');title('响应信号');024681012141618200.511.5冲击响应Time (sec)A m p l i t u d e024681012141618205101520阶跃响应Time (sec)A m p l i t u d e02468101214161820-1-0.500.51Time A激励信号2468101214161820-10123TimeA响应信号三． 1.>> clear; b=[1,3]; a=[1,3,2]; t=0:0.08:8; e=[exp(-3*t)]; sys=tf(b,a); lsim(sys,e,t);0.10.20.30.40.50.60.70.80.91Linear Simulation ResultsTime (sec)A m p l i t u d e2.>> clear; b=[1,3]; a=[1,2,2]; t=0:0.08:8; sys=tf(b,a); step(sys)0.20.40.60.811.21.41.6Step ResponseTime (sec)A m p l i t u d e3．>> clear; b=[1,3]; a=[1,2,1]; t=0:0.08:8; e=[exp(-2*t)]; sys=tf(b,a); lsim(sys,e,t);0.10.20.30.40.50.60.70.80.91Linear Simulation ResultsTime (sec)A m p l i t u d eDoc: 1.>> clear; B=[1]; A=[1,1,1]; sys=tf(B,A,-1); n=0:200;e=5+cos(0.2*pi*n)+2*sin(0.7*pi*n); r=lsim(sys,e); stem(n,r);0204060801001201401601802002.>> clear;B=[1,1,1];A=[1,-0.5,-0.5];sys=tf(B,A,-1);e=[1,zeros(1,100)];n=0:100;r=lsim(sys,e);stem(n,r);010203040506070809010021 / 21。
信号处理实验报告
一、实验目的本次实验旨在通过MATLAB软件平台，对数字信号处理的基本概念、原理和方法进行学习和实践。
通过实验，加深对以下内容的理解：1. 离散时间信号的基本概念和性质；2. 离散时间系统及其特性；3. 离散傅里叶变换（DFT）及其性质；4. 离散傅里叶逆变换（IDFT）及其应用；5. 窗函数及其在信号处理中的应用。
二、实验内容1. 离散时间信号的产生与性质（1）实验步骤：1.1 利用MATLAB生成以下离散时间信号：- 单位脉冲序列：δ[n]；- 单位阶跃序列：u[n]；- 矩形序列：R[n]；- 实指数序列：a^n；- 复指数序列：e^(jωn)。
1.2 分析并比较这些信号的性质，如自相关函数、功率谱密度等。
（2）实验结果：实验结果显示，不同类型的离散时间信号具有不同的性质。
例如，单位脉冲序列的自相关函数为δ[n]，功率谱密度为无穷大；单位阶跃序列的自相关函数为R[n]，功率谱密度为有限值；矩形序列的自相关函数为R[n]，功率谱密度为无穷大；实指数序列和复指数序列的自相关函数和功率谱密度均为有限值。
2. 离散时间系统及其特性（1）实验步骤：2.1 利用MATLAB构建以下离散时间系统：- 线性时不变系统：y[n] = x[n] a^n；- 非线性时不变系统：y[n] = x[n]^2；- 线性时变系统：y[n] = x[n] (1 + n)。
2.2 分析并比较这些系统的特性，如稳定性、因果性、线性时不变性等。
（2）实验结果：实验结果显示，不同类型的离散时间系统具有不同的特性。
例如，线性时不变系统的输出与输入之间存在线性关系，且满足时不变性；非线性时不变系统的输出与输入之间存在非线性关系，但满足时不变性；线性时变系统的输出与输入之间存在线性关系，但满足时变性。
3. 离散傅里叶变换（DFT）及其性质（1）实验步骤：3.1 利用MATLAB对以下离散时间信号进行DFT变换：- 单位脉冲序列：δ[n]；- 单位阶跃序列：u[n]；- 矩形序列：R[n]。
信号实验报告南理工
本次实验旨在通过实际操作加深对信号处理基本理论的理解，掌握信号频谱分析的方法，学习不同窗函数对信号频谱的影响，以及采样定理在信号处理中的应用。
通过实验，培养学生动手能力、分析问题和解决问题的能力。
二、实验原理1. 信号频谱分析：利用傅里叶变换将信号从时域转换为频域，分析信号的频率成分和能量分布。
2. 窗函数：在信号截取过程中，窗函数用于减少截取信号边缘的泄漏效应，提高频谱分析的准确性。
3. 采样定理：奈奎斯特采样定理指出，为了无失真地恢复原信号，采样频率应大于信号最高频率的两倍。
三、实验设备与软件1. 实验设备：示波器、信号发生器、计算机等。
2. 实验软件：MATLAB、Simulink等。
四、实验内容1. 信号频谱分析：（1）定义一个离散信号x[n]，计算其频谱X[k]。
（2）分别采用矩形窗、汉宁窗、汉明窗对信号进行截取，计算截取信号的频谱。
（3）比较不同窗函数对信号频谱的影响。
2. 采样定理验证：（1）根据奈奎斯特采样定理，确定信号的最大采样间隔和最小采样点数。
（2）通过改变采样点数，观察频谱变化，验证采样定理。
3. 周期性信号的DFT分析：（1）计算信号x[n]的周期T。
（2）通过补零和截取信号，分析周期性信号的DFT。
1. 在MATLAB中定义离散信号x[n]，并计算其频谱X[k]。
2. 分别采用矩形窗、汉宁窗、汉明窗对信号进行截取，计算截取信号的频谱。
3. 比较不同窗函数对信号频谱的影响。
4. 根据奈奎斯特采样定理，确定信号的最大采样间隔和最小采样点数。
5. 改变采样点数，观察频谱变化，验证采样定理。
6. 计算信号x[n]的周期T，通过补零和截取信号，分析周期性信号的DFT。
六、实验结果与分析1. 信号频谱分析：通过实验，发现不同窗函数对信号频谱的影响不同。
矩形窗频谱泄漏严重，汉宁窗和汉明窗能较好地抑制泄漏。
2. 采样定理验证：实验结果表明，当采样点数小于最小采样点数时，频谱发生严重混叠；当采样点数等于最小采样点数时，频谱能够无失真地恢复原信号。
东南大学统计信号处理实验实验三
统计信号处理实验三目的：掌握卡尔曼滤波滤波器的原理；内容：用雷达跟踪目标，目标的运动可以看成是在径向和横向内的二维运动，其运动方程和观测方程分别为：1111122222(1)()0100(1)()()0100(1)()0001(1)()()0001s k s k T v k v k u k s k s k T v k v k u k +⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥+⎣⎦⎣⎦⎣⎦⎣⎦11112222(1)(1)(1)(1)1000(1)(1)(1)0010(1)s k y k v k w k y k s k w k v k +⎡⎤⎢⎥+++⎡⎤⎡⎤⎡⎤⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥+++⎣⎦⎣⎦⎣⎦⎢⎥+⎣⎦1()s k 、1()v k 和1()y k 分别为径向距离、速度和观测值，而2()s k 、2()v k 和2()y k 分别为横向距离、速度和观测值。
1()u k 和2()u k 是状态噪声，是目标速度的波动；1()w k 和2()w k 是观测噪声；四种噪声的均值都为0，呈高斯分布，互不相关。
T 是雷达扫描一次的时间，此处设为1.0秒。
假设目标距离雷达约160Km 左右，径向初速度设为300 m/s ，并且在向雷达靠近，横向初速度设为0 m/s 。
这样它的径向速度波动大，而横向速度波动小，所以我们假设1()u k 的方差21u σ为300m/s ，2()u k 的方差22u σ为11.210-⨯m/s 。
鉴于雷达的观测误差，我们假设观测噪声1()w k 和2()w k 的方差21w σ和22w σ均为1.0Km 。
其中21u σ，22u σ，21w σ和22w σ的初始值不是最佳的，学生完全可自己修改以上参数，并观察计算结果的变化，给出最好的滤波效果。
任务：1） 试用αβ-滤波法对信号进行处理，并通过计算机模拟对其跟踪过程进行验证；2） 试求其Kalman 滤波方程，并通过计算机模拟对其跟踪过程进行验证；3） 假设目标在运动过程中发生了机动（速度在某个时刻突然发生了改变），试观测此时的αβ-滤波和Kalman 滤波结果，并对结果进行解释。
东南大学统计信号处理实验一
（1）似然比检测方法
Bayeｓ风险
0.88５5
0.21４0
0．11４5
0.5４24
0．８425
0.1581
0.１576
0．473８
０.７899
0.11６2
0.2101
0.4４24
检测概率 =n0／M；虚警概率 =n2/M；漏警概率 =ｎ1／M；
Ｂaｙes风险 = ＝
５、用相同的方法,通过改变判决的门限，观察检测方法的 、 、 和Ｂａyｅｓ风险的变化。
６、用相同的方法，通过改变噪声的方差,观察检测方法的 、 、 和Ｂayｅｓ风险的变化。
7、设计匹配滤波器h（t)=c*s(T-t)，通过使待检测信号x(ｔ)经过匹配滤波器，即和h（t）进行卷积,得到滤波以后的输出X(t）。
5、改变噪声的方差,观察检测方法的 、 、 和Bayｅs风险的变化;
6、将信号取样间隔减小一倍(相应的取样点数增加一倍)，观察似然比检测方法的 、 、 和Bａyes风险的变化；
7、根据 设计一个离散匹配滤波器，并观察 经过该滤波器以后的输出。
三、实验要求:
1、设计仿真计算的Matlab程序,给出软件清单;
采用似然比检测方法得到的仿真结果如下：
pｄ=0.8855,pf＝0.２1４０,pm=0．11４5，r=０.54２4。
利用基于Bayes准则的检测方法得到的仿真结果如下：
Ｐd=０.803２,Ｐf=０．１２64,Pm＝0.19６8,ｒ=０.449６。
比较可得：
采用似然比检测方法得到的检测概率较大,漏警概率较小；基于Bayes准则的检测方法得到的虚警概率较小，风险系数较小。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
统计信号处理实验三目的：掌握卡尔曼滤波滤波器的原理；内容：用雷达跟踪目标，目标的运动可以看成是在径向和横向内的二维运动，其运动方程和观测方程分别为：1111122222(1)()0100(1)()()0100(1)()0001(1)()()0001s k s k T v k v k u k s k s k T v k v k u k +⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥+⎣⎦⎣⎦⎣⎦⎣⎦11112222(1)(1)(1)(1)1000(1)(1)(1)0010(1)s k y k v k w k y k s k w k v k +⎡⎤⎢⎥+++⎡⎤⎡⎤⎡⎤⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥+++⎣⎦⎣⎦⎣⎦⎢⎥+⎣⎦1()s k 、1()v k 和1()y k 分别为径向距离、速度和观测值，而2()s k 、2()v k 和2()y k 分别为横向距离、速度和观测值。
1()u k 和2()u k 是状态噪声，是目标速度的波动；1()w k 和2()w k 是观测噪声；四种噪声的均值都为0，呈高斯分布，互不相关。
T 是雷达扫描一次的时间，此处设为1.0秒。
假设目标距离雷达约160Km 左右，径向初速度设为300 m/s ，并且在向雷达靠近，横向初速度设为0 m/s 。
这样它的径向速度波动大，而横向速度波动小，所以我们假设1()u k 的方差21u σ为300m/s ，2()u k 的方差22u σ为11.210-⨯m/s 。
鉴于雷达的观测误差，我们假设观测噪声1()w k 和2()w k 的方差21w σ和22w σ均为1.0Km 。
其中21u σ，22u σ，21w σ和22w σ的初始值不是最佳的，学生完全可自己修改以上参数，并观察计算结果的变化，给出最好的滤波效果。
任务：1） 试用αβ-滤波法对信号进行处理，并通过计算机模拟对其跟踪过程进行验证；2） 试求其Kalman 滤波方程，并通过计算机模拟对其跟踪过程进行验证；3） 假设目标在运动过程中发生了机动（速度在某个时刻突然发生了改变），试观测此时的αβ-滤波和Kalman 滤波结果，并对结果进行解释。
要求：1）设计仿真计算的Matlab 程序，给出软件清单；2）完成实验报告，给出实验结果，并对实验数据进行分析。
(1)αβ-滤波法对信号进行处理。
clear;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=[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)=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=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-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;取值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 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/(T^2) 0 0;0 0 sigma_w sigma_w/T;0 0 sigma_w/T sigma_u2+2*sigma_w/(T^2) ];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)=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==1P1=A*P2*A'+Q;elseP1=A*P0*A'+Q;endK=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('真实值','估计值');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)假设目标在运动过程中发生了机动（速度在某个时刻突然发生了改变），试观测此时的αβ-滤波和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=[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==199X(2,i+1)= X(2,i+1)+10;X(4,i+1)= X(4,i+1)+10;endY(:,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=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-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;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 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/(T^2) 0 0;0 0 sigma_w sigma_w/T;0 0 sigma_w/T sigma_u2+2*sigma_w/(T^2) ];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)];if i==199X(2,i+1)= X(2,i+1)+8;X(4,i+1)= X(4,i+1)+8;endY(:,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==1P1=A*P2*A'+Q;elseP1=A*P0*A'+Q;endK=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('真实值','估计值');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;。