实验报告-卡尔曼滤波.doc

数字信号处理实验报告
姓名： 专业: 通信与信息系统 学号: 日期：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，自编卡尔曼滤波递推程序，估计信号的波形。

任务二：

设计一维纳滤波器。
（1） 产生三组观测数据：首先根据产生信号，将其加噪（信噪比分别为20dB，10dB，6dB），得到观测数据，，。
（2） 估计，，，的AR模型参数。假设信号长度为，AR模型阶数为，分析实验结果，并讨论改变，对实验结果的影响。
实验任务一
1. 卡尔曼滤波原理
1.1 卡尔曼滤波简介
早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于空间技术的发展，为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推最优估计理论。它用状态空间法描述系统，由状态方程和量测方程所组成，即知道前一个状态的估计值和最近一个观测数据，采用递推的算法估计当前的状态值。由于卡尔曼滤波采用递推法，适合于计算机处理，并且可以用来处理多维和非平稳随机信号，现已广泛应用于很多领域，并取得了很好的结果。卡尔曼滤波一经出现，就受到人们的很大重视，并 在实践中不断丰富和完善，其中一个成功的应用是设计运载体的高精度组合导航系统。卡尔曼滤波具有以下的特点:
（1） 算法是递推的，且状态空间法采用在时域内设计滤波器的方法，因而适用于多维随机过程的估计；离散型卡尔曼算法适用于计算机处理。
（2） 用递推法计算，不需要知道全部过去的值，用状态方程描述状态变量的动态变化规律，因此信号可以是平稳的，也可以是非平稳的，即卡尔曼滤波适用于非平稳过程。
（3） 卡尔曼滤波采取的误差准则仍为估计误差的均方值最小。
1.2 卡尔曼滤波的状态方程和测量方程
假设某系统时刻的状态变量为，状态方程和量测方程（输出方程）表示为

其中，是状态变量；表示输入信号是白噪声；是观测噪声；是观测数据。
为了推导简单，假设状态变量的增益矩阵不随时间发生变化，，都是均值为零的正态白噪声，方差分别是和，并且初始状态与，都不相关，表示相关系数。即：

其中

1.3 卡尔曼滤波的递推算法
卡尔曼滤波采用递推算法来实现，其基本思想是先不考虑输入信号和观测噪声的影响，得到状态变量和输出信号（即观测数据）的估计值，再用输出信号的估计误差加权后校正状态变量的估计值，使状态变量估计误差的均方值最小。因此，卡尔曼滤波器的关键是计算出加权矩阵的最佳值。
当不考虑观测噪声和输入信号时，状态方程和量测方程为

显然，由于不考虑观测噪声的影响，输出信号的估计值与实际值是有误差的，用表示

为了提高状态估计的质量，用输出信号的估计误差来校正状态变量

其中，为增益矩阵，即加权矩阵。经过校正后的状态变量的估计误差及其均方值分别用和表示，把未经校正的状态变量的估计误差的均方值用表示

卡尔曼滤波要求状态变量的估计误差的均方值为最小，因此卡尔曼滤波的关键即为通过选择合适的，使得取得最小值。首先推导状态变量的估计值和状态变量的的估计误差，然后计算的均方值，通过化简，得到一组卡尔曼滤波的递推公式：

假设初始条件，，，，，，已知，其中，，那么递推流程如下：
，，
2. 卡尔曼滤波递推程序编程思想
题目分析
（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);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'*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);
end
end
M=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;
4. 实验结果

实验任务二
1. 维纳滤波器原理
维纳-霍夫方程

当是一个长度为的因果序列（即一个长度为的滤波器）时，维纳-霍夫方程表述为

定义

则可写成矩阵的形式，即

对上式求逆，得到

由以上式子可知：若已知期望信号与观测数据的互相关函数及观测数据的自相关函数，则可以通过矩阵求逆运算，得到维纳滤波器的最佳解。同时可以看到，直接从时域求解因果的维纳滤波器，当选择的滤波器的长度较大时，计算工作量很大，并且需要计算的逆矩阵，从而要求的存储量也很大

预测是根据观测到对的过去数据来估计当前或将来的信号值。维纳预测是已知以前时刻的个数据，估计当前时刻，或者未来时刻的信号值，即估计，估计得到的结果仍然要求满足均方误差最小的准则。信号可以预测是由于信号内部存在着关联性。预测是利用数据前后的关联性，根据其中一部分推知其余部分。
一步线性预测的时域解
已知，，…，，预测，假设噪声，这样的预测成为一步线性预测。设定系统的单位脉冲响应为。根据现行系统的基本理论，输出信号

令，则

预测误差

其中

要使均方误差为最小值，要求
，，...，
又因为，我们可以得到
，，...，
所以
，，...， （1）
由于预测器的输出是输入信号的线性组合，所以可得：

以上说明误差信号与输入信号满足正交性原理，预测误差与预测信号值同样满足正交性原理。预测误差的最小均方值
（2）
由（1）（2）联立方程组，写成矩阵形式可得

这就是有名的Yule-Walker（维纳-霍夫）方程。
2. 实验编程思想
在本实验中，首先根据要求产生加噪不同的观测数据，，，然后可利用已知条件代入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)); 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);
end
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)');ylabel('幅度');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 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);
sum = sum + X(j) * X(j + i - 1);
end
rx(i) = sum / (L - i + 1);
end
for i = 1 : N + 1
R(i, 1:(i-1)) = rx((i-1):-1:1);
R(i, i:(N+1)) = rx(1:(N - i + 2));
end
zx = 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)]);

function 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= awgn(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)');ylabel('幅度');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 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);
sum = sum + X(j) * X(j + i - 1);
end
rx(i) = sum / (L - i + 1);
end
for i = 1 : N + 1
R(i, 1:(i-1)) = rx((i-1):-1:1);
R(i, i:(N+1)) = rx(1:(N - i + 2));
end
zx = 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)]);

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信号来说 N 阶模型参数和误差分别为：
AR系数： 1 -0.26441
均方误差：1.0531
N=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.92331
N=3
对X1信号来说 N 阶模型参数和误差分别为：
AR系数： 1 -0.10673 0.050928 -0.19364
均方误差：1.4197
对X2信号来说 N 阶模型参数和误差分别为：
AR系数： 1 -0.35451 0.62013 -0.75585
均方误差：0.95739
对X3信号来说 N 阶模型参数和误差分别为：
AR系数： 1 -0.12221 0.14428 -0.34185
均方误差：0.99317
N=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.028944 0.21289 -0.2508
均方误差：1.3557
对X3信号来说 N 阶模型参数和误差分别为：
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.022914 0.0078698
均方误差：1.2033
对X2信号来说 N 阶模型参数和误差分别为：
AR系数： 1 0.017414 0.19629
均方误差：1.1607
对X3信号来说 N 阶模型参数和误差分别为：
AR系数： 1 -0.012789 0.13086
均方误差：1.1483
L=200
对X1信号来说 N 阶模型参数和误差分别为：
AR系数： 1 -0.17679 0.073726
均方误差：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=1

X1信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.15954
预测误差的最小均方值：1.0612
X2信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.1682
预测误差的最小均方值：1.1551
X3信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.1161
预测误差的最小均方值：1.2883
N=2
X1信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.20658 0.25733
预测误差的最小均方值：0.98824
X2信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.10574 0.15188
预测误差的最小均方值：1.0349
X3信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.17376 0.23089
预测误差的最小均方值：1.2323
N=5
X1信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.21587 0.22397 -0.24306 0.24469 -0.13453
预测误差的最小均方值：0.88869
X2信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.2537 0.31482 -0.19014 0.12243 0.040983
预测误差的最小均方值：0.89859
X3信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.27812 0.33384 -0.27881 0.25752 -0.11447
预测误差的最小均方值：1.0384
N=10
X1信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.085746 0.17042 -0.24887 0.3049 -0.29013 0.34871 -0.33129 0.48704 -0.39942 0.37265
预测误差的最小均方值：0.96799
X2信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.10183 0.20689 -0.18967 0.17418 -0.023448 -0.011321 0.0093937 0.18436 -0.1808 0.13523
预测误差的最小均方值：0.98478
X3信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.21104 0.51543 -0.60156 0.76905 -0.82715 0.92023 -0.90162 0.95174 -0.74381 0.6246
预测误差的最小均方值：0.96821

L=1000

X1信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.29706 0.2975 -0.25122 0.30972 -0.26171
预测误差的最小均方值：0.99801
X2信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.23764 0.2173 -0.15246 0.2091 -0.16053
预测误差的最小均方值：1.0381
X3信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.29651 0.28487 -0.21534 0.26745 -0.21993
预测误差的最小均方值：1.0984
L=500
X1信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.34838 0.37975 -0.33385 0.3446 -0.2554
预测误差的最小均方值：1.1584
X2信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.32197 0.34402 -0.24965 0.26156 -0.26372
预测误差的最小均方值：1.1777
X3信号: N 阶模型参数和预测误差的最小均方值分别为：
AR系数： 1 -0.29053 0.28967 -0.19091 0.17643 -0.14395
预测误差的最小均方值：1.405