实验报告-卡尔曼滤波

1、精选优质文档-倾情为你奉上数字信号处理实验报告姓名： 专业: 通信与信息系统 学号: 日期：2015.11实验内容任务一：一连续平稳的随机信号，自相关函数，信号为加性噪声所干扰，噪声是白噪声，测量值的离散值为已知，-3.2，-0.8，-14，-16，-17，-18，-3.3，-2.4，-18，-0.3，-0.4，-0.8，-19，-2.0，-1.2，-11，-14，-0.9，-0.8，10，0.2，0.5，-0.5，2.4，-0.5，0.5，-13，0.5，10，-12，0.5，-0.6，-15，-0.7，15，0.5，-0.7，-2.0，-19，-17，-11，-14，自编卡尔曼滤波递推程序

2、，估计信号的波形。任务二：设计一维纳滤波器。（1） 产生三组观测数据：首先根据产生信号，将其加噪（信噪比分别为20dB，10dB，6dB），得到观测数据，。（2） 估计，的AR模型参数。假设信号长度为，AR模型阶数为，分析实验结果，并讨论改变，对实验结果的影响。实验任务一1. 卡尔曼滤波原理1.1 卡尔曼滤波简介早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于空间技术的发展，为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推最优估计理论。它用状态空间法描述系统，由状态方程和量测方程所组成，即知道前一个状态的估计值和最近一个观测数据，采用递推的算法估计

3、当前的状态值。由于卡尔曼滤波采用递推法，适合于计算机处理，并且可以用来处理多维和非平稳随机信号，现已广泛应用于很多领域，并取得了很好的结果。卡尔曼滤波一经出现，就受到人们的很大重视，并 在实践中不断丰富和完善，其中一个成功的应用是设计运载体的高精度组合导航系统。卡尔曼滤波具有以下的特点:（1） 算法是递推的，且状态空间法采用在时域内设计滤波器的方法，因而适用于多维随机过程的估计；离散型卡尔曼算法适用于计算机处理。（2） 用递推法计算，不需要知道全部过去的值，用状态方程描述状态变量的动态变化规律，因此信号可以是平稳的，也可以是非平稳的，即卡尔曼滤波适用于非平稳过程。（3） 卡尔曼滤波采取的误差准

4、则仍为估计误差的均方值最小。1.2 卡尔曼滤波的状态方程和测量方程假设某系统时刻的状态变量为，状态方程和量测方程（输出方程）表示为其中，是状态变量；表示输入信号是白噪声；是观测噪声；是观测数据。为了推导简单，假设状态变量的增益矩阵不随时间发生变化，都是均值为零的正态白噪声，方差分别是和，并且初始状态与，都不相关，表示相关系数。即：其中1.3 卡尔曼滤波的递推算法卡尔曼滤波采用递推算法来实现，其基本思想是先不考虑输入信号和观测噪声的影响，得到状态变量和输出信号（即观测数据）的估计值，再用输出信号的估计误差加权后校正状态变量的估计值，使状态变量估计误差的均方值最小。因此，卡尔曼滤波器的关键是计算出

5、加权矩阵的最佳值。当不考虑观测噪声和输入信号时，状态方程和量测方程为显然，由于不考虑观测噪声的影响，输出信号的估计值与实际值是有误差的，用表示为了提高状态估计的质量，用输出信号的估计误差来校正状态变量其中，为增益矩阵，即加权矩阵。经过校正后的状态变量的估计误差及其均方值分别用和表示，把未经校正的状态变量的估计误差的均方值用表示卡尔曼滤波要求状态变量的估计误差的均方值为最小，因此卡尔曼滤波的关键即为通过选择合适的，使得取得最小值。首先推导状态变量的估计值和状态变量的的估计误差，然后计算的均方值，通过化简，得到一组卡尔曼滤波的递推公式：假设初始条件，已知，其中，那么递推流程如下：，2. 卡尔曼滤波

6、递推程序编程思想题目分析（1） 由于信号为加性噪声所干扰，可知，所以（2） 又因为噪声为白噪声，所以（3） 因为，所以由此可知，即，可得到：，因为抽样间隔，所以：。 （4）因此，所以因此编程分析由上面的分析可知初始条件，已知，在仿真中假设，则，由以上参数可得卡尔曼实际递推公式将得到的公式代入前面分析的递推公式，即可进行迭代得到结果。3. MATLAB源代码根据以上分析，编写matlab程序如下：%-卡尔曼滤波-%-说明%X(k+1)=Ak*X(k)+W(k);%Y(k)=Ck*X(k)+V(k)%clear;clc;%基本参数值Ak=exp(-0.02);Ck=1;Qk=1-exp(-0.04

7、);Rk=1;%初始值设置X0=0;P0=1;%观测值y(k)Y=-3.2 -0.8 -14 -16 -17 -18 -3.3 -2.4 -18 -0.3 -0.4 -0.8 -19 -2.0 -1.2 . -11 -14 -0.9 0.8 10 0.2 0.5 2.4 -0.5 0.5 -13 0.5 10 -12 0.5 -0.6 -15 -0.7 15 . 0.5 -0.7 -2.0 -19 -17 -11 -14;%数据长度N=length(Y);for k=1:N if k=1 %k=1时由初值开始计算 P_(k)=Ak*P0*Ak'+Qk; H(k)=P_(k)*Ck

8、9;*inv(Ck*P_(k)*Ck'+Rk); X(k)=Ak*X0+H(k)*(Y(k)-Ck*Ak*X0); I=eye(size(H(k); P(k)=(I-H(k)*Ck)*P_(k); else %k>1时，开始递推 %递推公式 P_(k)=Ak*P(k-1)*Ak'+Qk; H(k)=P_(k)*Ck'*inv(Ck*P_(k)*Ck'+Rk); X(k)=Ak*X(k-1)+H(k)*(Y(k)-Ck*Ak*X(k-1); I=eye(size(H(k); P(k)=(I-H(k)*Ck)*P_(k); endendM=1:N;T=0.02

9、*M;%作图，画出x(t)的波形figure(1)plot(T,Y,'r','LineWidth',1);xlabel('t');ylabel('y(t)');title('卡尔曼滤波-测量信号y(t)波形');grid;figure(2)plot(T,X,'b','LineWidth',1);xlabel('t');ylabel('x(t)');title('卡尔曼滤波-估计信号x(t)波形');grid;4. 实验结果实验任务二1.

10、 维纳滤波器原理维纳-霍夫方程当是一个长度为的因果序列（即一个长度为的滤波器）时，维纳-霍夫方程表述为定义则可写成矩阵的形式，即对上式求逆，得到由以上式子可知：若已知期望信号与观测数据的互相关函数及观测数据的自相关函数，则可以通过矩阵求逆运算，得到维纳滤波器的最佳解。同时可以看到，直接从时域求解因果的维纳滤波器，当选择的滤波器的长度较大时，计算工作量很大，并且需要计算的逆矩阵，从而要求的存储量也很大预测是根据观测到对的过去数据来估计当前或将来的信号值。维纳预测是已知以前时刻的个数据，估计当前时刻，或者未来时刻的信号值，即估计，估计得到的结果仍然要求满足均方误差最小的准则。信号可以预测是由于信号

11、内部存在着关联性。预测是利用数据前后的关联性，根据其中一部分推知其余部分。一步线性预测的时域解已知，预测，假设噪声，这样的预测成为一步线性预测。设定系统的单位脉冲响应为。根据现行系统的基本理论，输出信号令，则预测误差其中要使均方误差为最小值，要求，.，又因为，我们可以得到，.，所以，.， （1）由于预测器的输出是输入信号的线性组合，所以可得：以上说明误差信号与输入信号满足正交性原理，预测误差与预测信号值同样满足正交性原理。预测误差的最小均方值 （2）由（1）（2）联立方程组，写成矩阵形式可得这就是有名的Yule-Walker（维纳-霍夫）方程。2. 实验编程思想在本实验中，首先根据要求产生加噪

12、不同的观测数据，然后可利用已知条件代入Yule-Walker方程，即可求解AR模型参数。在本实验中，假设，信号的初值。3. MATLAB代码function Wiener_predict(L,N)%clc;clear;%信噪比SN1 = 6; SN2 = 10;SN3 = 20;%产生信号s(n)a=0.2;W = random('norm', 0, 1, L, 1); S(1) = 0;for n = 2 : L S(n) = a * S(n-1) + W(n); end%产生观测信号Am = sum(abs(S).2) / L;P1 = Am / (10(SN1/20);

13、P2 = Am / (10(SN2/20);P3 = Am / (10(SN3/20);V1 = random('norm', 0, P1, L, 1);V2 = random('norm', 0, P2, L, 1);V3 = random('norm', 0, P3, L, 1);for n=1:L X1(n) = S(n) + V1(n); X2(n) = S(n) + V2(n); X3(n) = S(n) + V3(n);endsubplot(2,2,1);plot(S,'b');title('信号S(n)

14、9;);ylabel('幅度');grid on;subplot(2,2,2);plot(X1,'b');title('观测信号X1(n)');ylabel('幅度');grid on;subplot(2,2,3);plot(X2,'b');title('观测信号X2(n)');ylabel('幅度');grid on;subplot(2,2,4);plot(X3,'b');title('观测信号X3(n)');ylabel('幅度')

15、;grid on;fprintf('n对X1信号来说 N 阶模型参数和误差分别为：n')AR(X1,N); fprintf('n对X2信号来说 N 阶模型参数和误差分别为：n')AR(X2,N); fprintf('n对X3信号来说 N 阶模型参数和误差分别为：n')AR(X3,N); function AR(X,N) L = length(X);rx = zeros(1, N + 1);R = zeros(N + 1, N + 1);for i = 1 : (N + 1) sum = 0; for j = 1 : (L - i + 1); su

16、m = sum + X(j) * X(j + i - 1); end rx(i) = sum / (L - i + 1);endfor i = 1 : N + 1 R(i, 1:(i-1) = rx(i-1):-1:1); R(i, i:(N+1) = rx(1:(N - i + 2);endzx = rx(2:(N+1);ap = inv(R(1:N, 1:N) * (-zx)'a = 1, ap'e = rx(1) + zx * ap;disp('AR系数： ',num2str(a);disp('均方误差：',num2str(e);funct

17、ion Wiener_new1(L,N)% 产生三组观测数据%信噪比(dB)SNR1 = 20; SNR2 = 10;SNR3 = 6;%产生信号s(n)a=0.2;W = random('norm', 0, 1, L, 1); S(1) = 0;for n = 2 : L S(n) = a * S(n-1) + W(n); end%加噪声产生观测限号X1= awgn(S,SNR1,'measured','linear'); X2= awgn(S,SNR2,'measured','linear'); X3= awg

18、n(S,SNR3,'measured','linear'); %画出信号图像subplot(2,2,1);plot(S,'b'); title('信号S(n)'); ylabel('幅度');grid on;subplot(2,2,2);plot(X1,'b');title('观测信号X1(n)');ylabel('幅度');grid on;subplot(2,2,3);plot(X2,'b');title('观测信号X2(n)');y

19、label('幅度');grid on;subplot(2,2,4);plot(X3,'b');title('观测信号X3(n)');ylabel('幅度');grid on;% 估计模型参数fprintf('n对X1信号来说 N 阶模型参数和误差分别为：n');AR(X1,N); fprintf('n对X2信号来说 N 阶模型参数和误差分别为：n');AR(X2,N); fprintf('n对X3信号来说 N 阶模型参数和误差分别为：n');AR(X3,N); function A

20、R(X,N) L = length(X);rx = zeros(1, N + 1);R = zeros(N + 1, N + 1);for i = 1 : (N + 1) sum = 0; for j = 1 : (L - i + 1); sum = sum + X(j) * X(j + i - 1); end rx(i) = sum / (L - i + 1);endfor i = 1 : N + 1 R(i, 1:(i-1) = rx(i-1):-1:1); R(i, i:(N+1) = rx(1:(N - i + 2);endzx = rx(2:(N+1);ap = inv(R(1:N,

21、 1:N) * (-zx)'a = 1, ap'e = rx(1) + zx * ap;disp('AR系数： ',num2str(a);disp('均方误差：',num2str(e);4. 实验结果与分析（1） 观测数据产生图1. 原始信号与观测信号（L=50）（2） 模型阶数N对实验结果的影响N=1对X1信号来说 N 阶模型参数和误差分别为：对X1信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.27766均方误差：1.1289对X2信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.29326均方误差：0.97283对X3

22、信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.26441均方误差：1.0531N=2对X1信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.34494 0.2854均方误差：1.0958对X2信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.1696 0.10742均方误差：1.1639对X3信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.19532 0.17033均方误差：0.92331N=3对X1信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.10673 0. -0.19364均方误差：1.4197对X2信号来说 N 阶模型参数和误

23、差分别为：AR系数： 1 -0.35451 0.62013 -0.75585均方误差：0.95739对X3信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.12221 0.14428 -0.34185均方误差：0.99317N=5对X1信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.35515 0.56619 -0.54005 0.65254 -0.51327均方误差：1.2405对X2信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.27343 0.10227 0. 0.21289 -0.2508均方误差：1.3557对X3信号来说 N 阶模型参数和误差分别为

24、：AR系数： 1 -0.36594 0.41414 -0.41665 0.66894 -0.60712均方误差：1.1025分析：由以上实验结果可知：在数据的长度一定的条件下，改变AR模型的阶数，均方误差会改变，当阶数在某个值时，均方误差的值最小，因此滤波器的阶数对实验结果有很大影响。在本次实验中，仿真情况有限，在以上仿真中我们可以看到当模型阶数N为某一固定值时，均方误差明显较小。（3） 信号长度L对实验结果的影响L=100对X1信号来说 N 阶模型参数和误差分别为：AR系数： 1 0. 0.均方误差：1.2033对X2信号来说 N 阶模型参数和误差分别为：AR系数： 1 0. 0.19629

25、均方误差：1.1607对X3信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0. 0.13086均方误差：1.1483L=200对X1信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.17679 0.均方误差：1.3371对X2信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.2649 0.16751均方误差：0.99844对X3信号来说 N 阶模型参数和误差分别为：AR系数： 1 -0.27145 0.17666均方误差：0.99289分析：由以上仿真结果可知，实验中存在误差，但仍然可以看出，随着信号长度的增加，均方误差减小，预测更准确。L=100，N=1X1信号

26、: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.15954预测误差的最小均方值：1.0612X2信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.1682预测误差的最小均方值：1.1551X3信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.1161预测误差的最小均方值：1.2883N=2X1信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.20658 0.25733预测误差的最小均方值：0.98824X2信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.10574 0.151

27、88预测误差的最小均方值：1.0349X3信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.17376 0.23089预测误差的最小均方值：1.2323N=5X1信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.21587 0.22397 -0.24306 0.24469 -0.13453预测误差的最小均方值：0.88869X2信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.2537 0.31482 -0.19014 0.12243 0.预测误差的最小均方值：0.89859X3信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0.27812 0.33384 -0.27881 0.25752 -0.11447预测误差的最小均方值：1.0384N=10X1信号: N 阶模型参数和预测误差的最小均方值分别为：AR系数： 1 -0. 0.17042 -0.24887 0.3