完整的实验报告下载连接
一、实验原理
卡尔曼滤波和维纳滤波都是最小均方误差为准则的线性估计器。卡尔曼滤波和维纳滤波的不同点在于：（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('平均误差');
