完整的实验报告下载连接https://download.csdn.net/download/LIsaWinLee/14884356

一、实验原理

    卡尔曼滤波和维纳滤波都是最小均方误差为准则的线性估计器。卡尔曼滤波和维纳滤波的不同点在于：（1）解决最佳滤波的方法不同，维纳滤波是用频域及传递函数的方法，卡尔曼是用时域及状态变量的方法；（2）维纳滤波要求过程的自相关系数和互相关函数的简单知识，而卡尔曼滤波则要求时域中状态变量及信号产生过程的详细知识；（3）维纳滤波要求平稳，而卡尔曼滤波则不要求。

1.1维纳滤波器原理

维纳滤波的原理是根据全部过去的和当前的观测数据xn,xn-1,… 来估计信号的当前值。以均方误差最小条件下求解系统的传递函数H(z) 或单位脉冲响应h(n) 。维纳滤波的信号模型是从信号与噪声的相关函数得到。

维纳滤波器的一般结构如下：

有一期望响应d(n) ，滤波器系数的设计准则是使得滤波器的输出y(n) 是均方意义上对期望响应的最优线性关系。

线性系统输出为：

yn=sn=mhmx(n-m)

均方误差为：Ee2n=Esn-m=0∞hmx(n-m)2

维纳滤波器的设计，实际上就是在最小均方误差的条件下，即∂E[e2n]∂hj=0 ，确定滤波器的冲激响应h(n) 或系统函数H(z) ，可等效于求解维纳-霍夫方程。

1.2卡尔曼滤波器原理

卡尔曼滤波的原理是不需要全部过去的观察数据，只根据前一个估计值和最近一个观察数据来估计信号的当前值。它是用状态空间法描述系统，即由状态方程和量测方程组成。解是以估计值的形式给出的，其算法是递推。卡尔曼滤波的信号模型（一维）如下：

离散系统的n维状态方程：x(k)=Akxk-1+wk-1

离散系统的m维测量方差：yk=ckxk+vk

令Ak 表示状态变量之间的增益矩阵，ck 为状态变量与输出信号之间的增益矩阵，不随时间发生变化，动态噪声wk 与观测噪声vk 都是零均值的正态噪声，且两者互不相关，Rk=var[vk] =E[vkvkT] 为量测噪声协方差矩阵，Qk=var[wk] =E[wkwkT] 为动态噪声协方差矩阵。

系统初始条件为Ex0=μ0

varx0 =E[( x-μ0)(x-μ0)T ]=p0

cov[x0 ,wk]=Ex0wkT=0

cov[x0 ,vk]=Ex0vkT=0

卡尔曼滤波的基本思想是先不考虑激励噪声wk 和观测噪声vk ，得到的状态估计值xk 和观测数据的估计值yk ，再用观测数据的估计误差yk=yk-yk 去修正状态的估计值xk ，通过选择修正矩阵H使得状态估计误差的均方值Pk 最小。

卡尔曼滤波的递推公式如下：

xk=Akxk-1+Hk(yk-CkAkxk-1)

Hk=Pk'CkT(CkPk'CkT+Rk)-1

Pk'=AkPk-1AkT+Qk-1

Pk=(I-HkCk)Pk'

假设初始条件Ak,Ck,Qk,Rk,yk ,xk-1,Pk-1 已知，其中x0=Ex0 , P0=var[x0] ，则卡尔曼滤波的递推流程如下

二、实验内容

随机信号 服从 过程，它是一个宽带过程，参数如下：

（1）通过观测方程是方差为1的高斯白噪声，要求分别利用Wiener滤波器和Kalman滤波器通过测量信号估计的波形。

（2）将 的方差改为4，按（1）中要求，重新估计 的波形。

自行选择Wiener滤波器的阶数，自行选择实验扩展内容，写出实验报告，并进行相关分析。

三、实验过程

3.1维纳滤波器

利用维纳滤波测量信号时，建立模型如下：

       1）将随机信号X(n)看成是由典型白噪声序列源W(n)激励一个线性系统产生，用一个差分方程表示xn-1.352xn-1+1.338xn-2-0.662xn-3+0.24xn-4=w(n) 。进行Z变换得到Hz=X(z)W(z)=11-1.352z-1+1.338z-2-0.662z-3+0.24z-4 ，均值为1的高斯白噪声序列W(n)可以用randn函数产生，再利用函数X=filter(B,A,W)产生随机信号X(n)。

2）观测方程是Y(n)=X(n)+V(n)，V(n)是方差为1的高斯白噪声，产生进入维纳滤波器的信号。

维纳滤波仿真流程如下：

3.2卡尔曼滤波器

卡尔曼滤波实际由两个过程组成：预测与校正。在预测阶段，滤波器使用上一状态的估计，做出对当前状态的预测。在校正阶段，滤波器利用对当前状态的观测值修正在预测阶段获得的预测值，以获得一个更接进真实值的新估计值。采用最陡下降算法，搜索方向为梯度负方向，每一步更新都使梯度函数值最小。

卡尔曼滤波仿真流程如下：

四、实验结果

4.1 维纳滤波器实验结果

（1）维纳滤波器的阶数取101阶，观测点数取100，噪声方差为1时，通过仿真得到真实轨迹，观测样本及估计轨迹的比较图像，和通过仿真得到的平均误差如图4.1所示。

        (a)真实轨迹、观测样本及估计轨迹的比较                                                            (b)平均误差

                                                                  图4.1 维纳滤波（噪声方差取1）

从图4.1(a)可以看出，利用维纳滤波器通过测量信号估计真实信号的波形，得到的估计轨迹接近于真实轨迹。当噪声方差为1时，估计轨迹与真实轨迹间的误差很小。

从图4.1(b)可以看出，利用维纳滤波器通过测量信号估计真实信号的波形，当噪声方差为1时，得到的估计轨迹与真实轨迹的平均误差较小，而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐减小的趋势。

（2）维纳滤波器的阶数取101阶，观测点数取100，噪声方差为4时，通过仿真得到真实轨迹，观测样本及估计轨迹的比较图像，和通过仿真得到的平均误差如图4.2所示。

               (a)真实轨迹、观测样本及估计轨迹的比较                                                          (b)平均误差

                                                                               图4.2 维纳滤波（噪声方差取4）

从图4.2(a)可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的误差变大了。

从图4.2(b)可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的平均误差变大了。而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐见效的趋势。

4.2 卡尔曼滤波器实验结果

（1）卡尔曼滤波器的阶数取101阶，观测点数取100，噪声方差为1时，通过仿真得到真实轨迹，观测样本及估计轨迹的比较图像，和通过仿真得到的平均误差如图4.1所示。

   

                   (a)真实轨迹、观测样本及估计轨迹的比较                                                (b)平均误差

                                                                    图4.3 卡尔曼滤波（噪声方差取1）

从图4.3(a)可以看出，利用卡尔曼滤波器通过测量信号估计真实信号的波形，得到的估计轨迹接近于真实轨迹。当噪声方差为1时，估计轨迹与真实轨迹间的误差很小。

从图4.3(b)可以看出，利用卡尔曼滤波器通过测量信号估计真实信号的波形，当噪声方差为1时，得到的估计轨迹与真实轨迹的平均误差较小，而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐减小的趋势。

（2）卡尔曼滤波器的阶数取101阶，观测点数取100，噪声方差为4时，通过仿真得到真实轨迹，观测样本及估计轨迹的比较图像，和通过仿真得到的平均误差如图4.4所示。

   

                    (a)真实轨迹、观测样本及估计轨迹的比较                                                             (b)平均误差

                                                                               图4.4 卡尔曼滤波（噪声方差取4）

从图4.4(a)可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的误差变大了。

从图4.4(b)可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的平均误差变大了。而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐见效的趋势。

五、总结与展望

通过此次实验，学会运用MATLAB工具实现构建维纳滤波器和卡尔曼滤波器，并对AR模型进行实验。从实验中，我了解到维纳滤波器是以LMS算法为核心的，而卡尔曼滤波器是从RLS算法演变而来的。从这两点则可以衍生出很多解题思路。最后通过MATLAB演示了信号在加入噪声后，通过维纳滤波器和卡尔曼滤波器进行估计，比较估计值和真实值，发现差别非常的小。希望后面若有遇到维纳滤波和卡尔曼滤波问时，能够回忆起课堂上和实验中学到的知识，充分的利用起来。

六、附录

6.1 维纳滤波器程序

（1）噪声方差为1时：

clc;

clear all;

maxlag=100;

N=100;%观测点数取100

x=zeros(N,1);

y=zeros(N,1);

%%列出状态方程%%

x(1)=randn(1,1);

x(2)=randn(1,1)+1.352*x(1);

x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);

x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);

for n=5:N

    x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);%x为真实值

end

v=randn(N,1);

y=x+v;%z_x为观测样本值=真值+噪声

%%滤波%%

x=x';

y=y';

xk_s(1)=y(1);%赋初值

xk_s(2)=y(2);

xk_s(3)=y(3);

xk_s(4)=y(4);

xk=[y(1);y(2);y(3);y(4)];

%%维纳滤波器的生成

[rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数

rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵，滤波器的阶数为101阶

rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数

rx2=rx2(101:end);

h=inv(rx1)*rx2';%维纳霍夫方程

xk_s=filter(h,1,y);%加噪信号通过滤波器后的输出

e_x=0;

eq_x=0;

e_x1=N:1;

%计算滤波的均值，计算滤波误差的均值

for i=1:N

    e_x(i)=x(i)-xk_s(i);%误差=真实值-滤波估计值

end

%%作图%%

t=1:N;

figure(1);

plot(t,x,'r-',t,y,'g:',t,xk_s,'b-.');

legend('真是轨迹','观测样本','估计轨迹');

figure(2);

plot(e_x);

legend('平均误差');

（2）噪声方差为4时：

clc;

clear all;

maxlag=100;

N=100;%观测点数取100

x=zeros(N,1);

y=zeros(N,1);

%%列出状态方程%%

x(1)=randn(1,1);

x(2)=randn(1,1)+1.352*x(1);

x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);

x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);

for n=5:N

    x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);%x为真实值

end

v=4*randn(N,1);

y=x+v;%z_x为观测样本值=真值+噪声

%%滤波%%

x=x';

y=y';

xk_s(1)=y(1);%赋初值

xk_s(2)=y(2);

xk_s(3)=y(3);

xk_s(4)=y(4);

xk=[y(1);y(2);y(3);y(4)];

%%维纳滤波器的生成

[rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数

rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵，滤波器的阶数为101阶

rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数

rx2=rx2(101:end);

h=inv(rx1)*rx2';%维纳霍夫方程

xk_s=filter(h,1,y);%加噪信号通过滤波器后的输出

e_x=0;

eq_x=0;

e_x1=N:1;

%计算滤波的均值，计算滤波误差的均值

for i=1:N

    e_x(i)=x(i)-xk_s(i);%误差=真实值-滤波估计值

end

%%作图%%

t=1:N;

figure(1);

plot(t,x,'r-',t,y,'g:',t,xk_s,'b-.');

legend('真是轨迹','观测样本','估计轨迹');

figure(2);

plot(e_x);

legend('平均误差');

6.2 卡尔曼滤波器程序

（1）噪声方差为1时：

N=100;

x=zeros(N,1);

y=zeros(N,1);

var=1;

I=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];

x(1)=randn(1,1);

x(2)=randn(1,1)+1.352*x(1);

x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);

x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);

for n=5:N

    x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);

end

v=randn(N,1);

y=x+v;

xk_s(1)=y(1);

xk_s(2)=y(2);

xk_s(3)=y(3);

xk_s(4)=y(4);

Ak=[1.352,-1.338,0.662,0.240;1,0,0,0;0,1,0,0;0,0,1,0];

Ck=[1,0,0,0];

Rk=[1];

Pk=[1 0 0 0

    0 1 0 0

    0 0 1 0

    0 0 0 1];

xk=[y(1);y(2);y(3);y(4)];

Qk=[1];

for r=5:N

    yk=y(r);

    Pk1=Ak*Pk*Ak'+Qk;

    Hk=Pk1*Ck'*inv(Ck*Pk1*Ck'+Rk);

    xk=Ak*xk+Hk*(yk-Ck*Ak*xk);

    Pk=(I-Hk*Ck)*Pk1;

    xk_s(r)=xk(1,1);

end

e_x=0;

for i=1:N

    e_x(i)=x(i)-xk_s(i);

end

t=1:N;

figure(1);

plot(t,x,'r-',t,y,'g:',t,xk_s,'b-');

legend('真实轨迹','观测样本','估计轨迹');

figure(2);

plot(e_x);

legend('平均误差');

（2）噪声方差为4时：

N=100;

x=zeros(N,1);

y=zeros(N,1);

var=2;

I=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];

x(1)=randn(1,1);

x(2)=randn(1,1)+1.352*x(1);

x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);

x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);

for n=5:N

    x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);

end

v=4*randn(N,1);

y=x+v;

xk_s(1)=y(1);

xk_s(2)=y(2);

xk_s(3)=y(3);

xk_s(4)=y(4);

Ak=[1.352,-1.338,0.662,0.240;1,0,0,0;0,1,0,0;0,0,1,0];

Ck=[1,0,0,0];

Rk=[1];

Pk=[1 0 0 0

    0 1 0 0

    0 0 1 0

    0 0 0 1];

xk=[y(1);y(2);y(3);y(4)];

Qk=[1];

for r=5:N

    yk=y(r);

    Pk1=Ak*Pk*Ak'+Qk;

    Hk=Pk1*Ck'*inv(Ck*Pk1*Ck'+Rk);

    xk=Ak*xk+Hk*(yk-Ck*Ak*xk);

    Pk=(I-Hk*Ck)*Pk1;

    xk_s(r)=xk(1,1);

end

e_x=0;

for i=1:N

    e_x(i)=x(i)-xk_s(i);

end

t=1:N;

figure(1);

plot(t,x,'r-',t,y,'g:',t,xk_s,'b-');

legend('真实轨迹','观测样本','估计轨迹');

figure(2);

plot(e_x);

legend('平均误差');