随机数字信号处理实验报告一——维纳滤波和卡尔曼滤波

完整的实验报告下载连接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('平均误差');

 

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值