实验报告卡尔曼滤波

数字信号解决实验报告姓名：专业:通信与信息系统学号:日期：.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，自编卡尔曼滤波递推程序，预计信号的波形。任务二：设计一维纳滤波器。产生三组观察数据：首先根据产生信号，将其加噪（信噪比分别为20dB，10dB，6dB），得到观察数据，，。预计，，，的AR模型参数。假设信号长度为，AR模型阶数为，分析实验成果，并讨论变化，对实验成果的影响。实验任务一卡尔曼滤波原理1.1卡尔曼滤波介绍早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于空间技术的发展，为理解决对非平稳、多输入输出随机序列的预计问题，卡尔曼提出了递推最优预计理论。它用状态空间法描述系统，由状态方程和量测方程所构成，即懂得前一种状态的预计值和近来一种观察数据，采用递推的算法预计现在的状态值。由于卡尔曼滤波采用递推法，适合于计算机解决，并且能够用来解决多维和非平稳随机信号，现已广泛应用于诸多领域，并获得了较好的成果。卡尔曼滤波一经出现，就受到人们的很大重视，并在实践中不停丰富和完善，其中一种成功的应用是设计运载体的高精度组合导航系统。卡尔曼滤波含有下列的特点:算法是递推的，且状态空间法采用在时域内设计滤波器的办法，因而合用于多维随机过程的预计；离散型卡尔曼算法合用于计算机解决。用递推法计算，不需要懂得全部过去的值，用状态方程描述状态变量的动态变化规律，因此信号能够是平稳的，也能够是非平稳的，即卡尔曼滤波合用于非平稳过程。卡尔曼滤波采用的误差准则仍为预计误差的均方值最小。1.2卡尔曼滤波的状态方程和测量方程假设某系统时刻的状态变量为，状态方程和量测方程（输出方程）表达为其中，是状态变量；表达输入信号是白噪声；是观察噪声；是观察数据。为了推导简朴，假设状态变量的增益矩阵不随时间发生变化，，都是均值为零的正态白噪声，方差分别是和，并且初始状态与，都不有关，表达有关系数。即：其中1.3卡尔曼滤波的递推算法卡尔曼滤波采用递推算法来实现，其基本思想是先不考虑输入信号和观察噪声的影响，得到状态变量和输出信号（即观察数据）的预计值，再用输出信号的预计误差加权后校正状态变量的预计值，使状态变量预计误差的均方值最小。因此，卡尔曼滤波器的核心是计算出加权矩阵的最佳值。当不考虑观察噪声和输入信号时，状态方程和量测方程为显然，由于不考虑观察噪声的影响，输出信号的预计值与实际值是有误差的，用表达为了提高状态预计的质量，用输出信号的预计误差来校正状态变量其中，为增益矩阵，即加权矩阵。通过校正后的状态变量的预计误差及其均方值分别用和表达，把未经校正的状态变量的预计误差的均方值用表达卡尔曼滤波规定状态变量的预计误差的均方值为最小，因此卡尔曼滤波的核心即为通过选择适宜的，使得获得最小值。首先推导状态变量的预计值和状态变量的的预计误差，然后计算的均方值，通过化简，得到一组卡尔曼滤波的递推公式：假设初始条件，，，，，，已知，其中，，那么递推流程以下：，，卡尔曼滤波递推程序编程思想题目分析由于信号为加性噪声所干扰，可知，因此又由于噪声为白噪声，因此由于，因此由此可知，，即，可得到：，由于抽样间隔，因此：。（4）因此，因此因此编程分析由上面的分析可知初始条件，，，，已知，在仿真中假设，则，，由以上参数可得卡尔曼实际递推公式将得到的公式代入前面分析的递推公式，即可进行迭代得到成果。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);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.90.8100.20.52.4-0.50.5-130.510-120.5-0.6-15-0.715...0.5-0.7-2.0-19-17-11-14];%数据长度N=length(Y);fork=1:Nifk==1%k=1时由初值开始计算P_(k)=Ak*P0*Ak'+Qk;H(k)=P_(k)*Ck'*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*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;实验成果实验任务二维纳滤波器原理维纳-霍夫方程当是一种长度为的因果序列（即一种长度为的滤波器）时，维纳-霍夫方程表述为定义则可写成矩阵的形式，即对上式求逆，得到由以上式子可知：若已知盼望信号与观察数据的互有关函数及观察数据的自有关函数，则能够通过矩阵求逆运算，得到维纳滤波器的最佳解。同时能够看到，直接从时域求解因果的维纳滤波器，当选择的滤波器的长度较大时，计算工作量很大，并且需要计算的逆矩阵，从而规定的存储量也很大预测是根据观察到对的过去数据来预计现在或将来的信号值。维纳预测是已知以前时刻的个数据，预计现在时刻，或者将来时刻的信号值，即预计，预计得到的成果仍然规定满足均方误差最小的准则。信号能够预测是由于信号内部存在着关联性。预测是运用数据前后的关联性，根据其中一部分推知其它部分。一步线性预测的时域解已知，，…，，预测，假设噪声，这样的预测成为一步线性预测。设定系统的单位脉冲响应为。根据现行系统的基本理论，输出信号令，则预测误差其中要使均方误差为最小值，规定，，...，又由于，我们能够得到，，...，因此，，...，（1）由于预测器的输出是输入信号的线性组合，因此可得：以上阐明误差信号与输入信号满足正交性原理，预测误差与预测信号值同样满足正交性原理。预测误差的最小均方值（2）由（1）（2）联立方程组，写成矩阵形式可得这就是有名的Yule-Walker（维纳-霍夫）方程。实验编程思想在本实验中，首先根据规定产生加噪不同的观察数据，，，然后可运用已知条件代入Yule-Walker方程，即可求解AR模型参数。在本实验中，假设，信号的初值。MATLAB代码functionWiener_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;forn=2:LS(n)=a*S(n-1)+W(n);end%产生观察信号Am=sum(abs(S).^2)/L;P1=Am/(10^(SN1/20));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);forn=1:LX1(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)');ylabel('幅度');gridon;subplot(2,2,2);plot(X1,'b');title('观察信号X1(n)');ylabel('幅度');gridon;subplot(2,2,3);plot(X2,'b');title('观察信号X2(n)');ylabel('幅度');gridon;subplot(2,2,4);plot(X3,'b');title('观察信号X3(n)');ylabel('幅度');gridon;fprintf('\n对X1信号来说N阶模型参数和误差分别为：\n')AR(X1,N);fprintf('\n对X2信号来说N阶模型参数和误差分别为：\n')AR(X2,N);fprintf('\n对X3信号来说N阶模型参数和误差分别为：\n')AR(X3,N);functionAR(X,N)L=length(X);rx=zeros(1,N+1);R=zeros(N+1,N+1);fori=1:(N+1)sum=0;forj=1:(L-i+1);sum=sum+X(j)*X(j+i-1);endrx(i)=sum/(L-i+1);endfori=1:N+1R(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)]);functionWiener_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;forn=2:LS(n)=a*S(n-1)+W(n);end%加噪声产生观察限号X1=awgn(S,SNR1,'measured','linear');X2=awgn(S,SNR2,'measured','linear');X3=awgn(S,SNR3,'measured','linear');%画出信号图像subplot(2,2,1);plot(S,'b');title('信号S(n)');ylabel('幅度');gridon;subplot(2,2,2);plot(X1,'b');title('观察信号X1(n)');ylabel('幅度');gridon;subplot(2,2,3);plot(X2,'b');title('观察信号X2(n)');ylabel('幅度');gridon;subplot(2,2,4);plot(X3,'b');title('观察信号X3(n)');ylabel('幅度');gridon;%%预计ＡＲ模型参数fprintf('\n对X1信号来说N阶模型参数和误差分别为：\n');AR(X1,N);fprintf('\n对X2信号来说N阶模型参数和误差分别为：\n');AR(X2,N);fprintf('\n对X3信号来说N阶模型参数和误差分别为：\n');AR(X3,N);functionAR(X,N)L=length(X);rx=zeros(1,N+1);R=zeros(N+1,N+1);fori=1:(N+1)sum=0;forj=1:(L-i+1);sum=sum+X(j)*X(j+i-1);endrx(i)=sum/(L-i+1);endfori=1:N+1R(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)]);实验成果与分析观察数据产生图1.原始信号与观察信号（L=50）模型阶数N对实验成果的影响N=1对X1信号来说N阶模型参数和误差分别为：对X1信号来说N阶模型参数和误差分别为：AR系数：1-0.27766均方误差：1.1289对X2信号来说N阶模型参数和误差分别为：AR系数：1-0.29326均方误差：0.97283对X3信号来说N阶模型参数和误差分别为：AR系数：1-0.26441均方误差：1.0531N=2对X1信号来说N阶模型参数和误差分别为：AR系数：1-0.344940.2854均方误差：1.0958对X2信号来说N阶模型参数和误差分别为：AR系数：1-0.16960.10742均方误差：1.1639对X3信号来说N阶模型参数和误差分别为：AR系数：1-0.195320.17033均方误差：0.92331N=3对X1信号来说N阶模型参数和误差分别为：AR系数：1-0.106730.050928-0.19364均方误差：1.4197对X2信号来说N阶模型参数和误差分别为：AR系数：1-0.354510.6-0.75585均方误差：0.95739对X3信号来说N阶模型参数和误差分别为：AR系数：1-0.122210.14428-0.34185均方误差：0.99317N=5对X1信号来说N阶模型参数和误差分别为：AR系数：1-0.355150.56619-0.540050.65254-0.51327均方误差：1.2405对X2信号来说N阶模型参数和误差分别为：AR系数：1-0.273430.102270.0289440.21289-0.2508均方误差：1.3557对X3信号来说N阶模型参数和误差分别为：AR系数：1-0.365940.41414-0.416650.66894-0.60712均方误差：1.1025分析：由以上实验成果可知：在数据的长度一定的条件下，变化AR模型的阶数，均方误差会变化，当阶数在某个值时，均方误差的值最小，因此滤波器的阶数对实验成果有很大影响。在本次实验中，仿真状况有限，在以上仿真中我们能够看到当模型阶数N为某一固定值时，均方误差明显较小。信号长度L对实验成果的影响L=100对X1信号来说N阶模型参数和误差分别为：AR系数：10.0229140.0078698均方误差：1.2033对X2信号来说N阶模型参数和误差分别为：AR系数：10.0174140.19629均方误差：1.1607对X3信号来说N阶模型参数和误差分别为：AR系数：1-0.0127890.13086均方误差：1.1483L=200对X1信号来说N阶模型参数和误差分别为：AR系数：1-0.176790.073726均方误差：1.3371对X2信号来说N阶模型参数和误差分别为：AR系数：1-0.26490.16751均方误差：0.99844对X3信号来说N阶模型参数和误差分别为：AR系数：1-0.271450.17666均方误差：0.99289分析：由以上仿真成果可知，实验中存在误差，但仍然能够看出，随着信号长度的增加，均方误差减小，预测更精确。L=100，N=1X1信号: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.206580.25733预测误差的最小均方值：0.98824X2信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.105740.15188预测误差的最小均方值：1.0349X3信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.173760.23089预测误差的最小均方值：1.2323N=5X1信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.215870.22397-0.243060.24469-0.13453预测误差的最小均方值：0.88869X2信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.25370.31482-0.190140.122430.040983预测误差的最小均方值：0.89859X3信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.278120.33384-0.278810.25752-0.11447预测误差的最小均方值：1.0384N=10X1信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.0857460.17042-0.248870.3049-0.290130.34871-0.331290.48704-0.399420.37265预测误差的最小均方值：0.96799X2信号:N阶模型参数和预测误差的最小均方值分别为：AR系数：1-0.101830.20689-0.18967