维纳、卡尔曼滤波简介及MATLAB实现

1、现代数字信号处理 课程作业维纳、卡尔曼、RLS、LMS算法matlab实现维纳滤波从噪声中提取信号波形的各种估计方法中，维纳（Wiener）滤波是一种最基本的方法，适用于需要从噪声中分离出的有用信号是整个信号（波形），而不只是它的几个参量。设维纳滤波器的输入为含噪声的随机信号。期望输出与实际输出之间的差值为误差，对该误差求均方，即为均方误差。因此均方误差越小，噪声滤除效果就越好。为使均方误差最小，关键在于求冲激响应。如果能够满足维纳霍夫方程，就可使维纳滤波器达到最佳。维纳滤波器的优点是适应面较广，无论平稳随机过程是连续的还是离散的，是标量的还是向量的，都可应用。维纳滤波器的缺点是，要求得到半无

2、限时间区间内的全部观察数据的条件很难满足，同时它也不能用于噪声为非平稳的随机过程的情况，对于向量情况应用也不方便。因此，维纳滤波在实际问题中应用不多。下面是根据维纳滤波器给出的图像处理matlab实例，在下面实例中维纳滤波和均值滤波相比较，并且做了维纳复原、边缘提取、图像增强的实验：%*维纳滤波和均值滤波的比较*I=imread(lena.bmp);J=imnoise(I,gaussian,0,0.01);Mywiener2 = wiener2(J,3 3);Mean_temp = ones(3,3)/9;Mymean = imfilter(J,Mean_temp); figure(1);su

3、bplot(121),imshow(Mywiener2),title(维纳滤波器输出);subplot(122),imshow(uint8(Mymean),),title(均值滤波器的输出);%*维纳复原程序*figure(2);subplot(231),imshow(I),title(原始图像); LEN = 20;THETA =10;PSF = fspecial(motion,LEN,THETA); Blurred = imfilter(I,PSF,circular);subplot(232),imshow(Blurred),title(生成的运动的模糊的图像);noise = 0.1*r

4、andn(size(I);subplot(233),imshow(im2uint8(noise),title(随机噪声);BlurredNoisy=imadd(Blurred,im2uint8(noise);subplot(234),imshow(BlurredNoisy),title(添加了噪声的模糊图像); Move=deconvwnr(Blurred,PSF);subplot(235),imshow(Move),title(还原运动模糊的图像);nsr = sum(noise(:).2)/sum(im2double(I(:).2);wnr2 = deconvwnr(BlurredNois

5、y,PSF,nsr);subplot(236),imshow(wnr2),title(还原添加了噪声的图像);%*维纳滤波应用于边缘提取*N = wiener2(I,3,3);%选用不同的维纳窗在此修改M = I - N;My_Wedge = im2bw (M,5/256);%化二值图像BW1 = edge(I,prewitt);BW2 = edge(I,canny);BW3 = edge(I,zerocross);BW4 = edge(I,roberts); figure(3)subplot(2,4,3 4 7 8),imshow(My_Wedge),title(应用维纳滤波进行边沿提取);

6、subplot(241),imshow(BW1),title(prewitt);subplot(242),imshow(BW2),title(canny);subplot(245),imshow(BW3),title(zerocross);subplot(246),imshow(BW4),title(roberts);%*维纳滤波应用于图像增强*for i = 1 2 3 4 5 K = wiener2(I,5,5);end K = K + I;figure(4);subplot(121),imshow(I),title(原始图像);subplot(122),imshow(K),title(增

7、强后的图像);卡尔曼滤波卡尔曼滤波的一个典型实例是从一组有限的，对物体位置的，包含噪声的观察序列预测出物体的坐标位置及速度。卡尔曼滤波利用目标的动态信息,设法去掉噪声的影响，得到一个关于目标位置的好的估计。这个估计可以是对当前目标位置的估计(滤波)，也可以是对于将来位置的估计(预测)，也可以是对过去位置的估计(插值或平滑)。状态估计是卡尔曼滤波的重要组成部分。一般来说，根据观测数据对随机量进行定量推断就是估计问题，特别是对动态行为的状态估计，它能实现实时运行状态的估计和预测功能。受噪声干扰的状态量是个随机量，不可能测得精确值，但可对它进行一系列观测，并依据一组观测值，按某种统计观点对它进行估计

8、。使估计值尽可能准确地接近真实值，这就是最优估计。真实值与估计值之差称为估计误差。卡尔曼提出的递推最优估计理论，采用状态空间描述法，在算法采用递推形式，卡尔曼滤波能处理多维和非平稳的随机过程。卡尔曼滤波的matlab实现如下所示：10%卡尔曼滤波器的应用，运动轨迹估计clc;clear;T=4;%雷达扫描周期num=100;%滤波次数%产生真实轨迹N=800/T;x=zeros(N,1);y=zeros(N,1);vx=zeros(N,1);vy=zeros(N,1);x(1)=-2000;y(1)=500;vx=10;vy=0;ax=0;ay=0;var=100;for i=1:N-1x(i

9、+1)=x(i)+vx*T+0.5*ax*T2;y(i+1)=y(i)+vy*T+0.5*ay*T2;endnx=zeros(N,1);ny=zeros(N,1);nx=100*randn(N,1);ny=100*randn(N,1);zx=x+nx;zy=y+ny;%滤波100次for m=1:numz=2:1;xks(1)=zx(1);yks(1)=zy(1);xks(2)=zx(2);yks(2)=zy(2);o=4:4;g=4:2;h=2:4; q=2:2;xk=4:1;perr=4:4;o=1,T,0,0;0,1,0,0;0,0,1,T;0,0,0,1;h=1 0 0 0;0 0 1

10、 0;g=T/2,0;T/2,0;0,T/2;0,T/2;q=10000 0;0 10000;perr=var2 var2/T 0 0var*var/T 2*var2/(T2) 0 00 0 var2 var2/T0 0 var2/T 2*var2/(T2);vx=(zx(2)-zx(1)/2;vy=(zy(2)-zy(1)/2;xk=zx(1);vx;zy(1);vy;%Kalman 滤波开始for r=3:N;z=zx(r);zy(r);xk1=o*xk;perr1=o*perr*o;k=perr1*h*inv(h*perr1*h+q);xk=xk1+k*(z-h*xk1);perr=(e

11、ye(4)-k*h)*perr1;xks(r)=xk(1,1);yks(r)=xk(3,1);vkxs(r)=xk(2,1);vkys(r)=xk(4,1);xk1s(r)=xk1(1,1);ykls(r)=xk1(3,1);perr11(r)=perr(1,1);perr12(r)=perr(1,2);perr22(r)=perr(2,2);rex(m,r)=xks(r);rey(m,r)=yks(r);end %结束一次滤波endex=0;ey=0;eqx=0;eqy=0;ey1=0;ex1=N:1;ey1=N:1;%计算滤波的均值，计算滤波误差的均值for i=1:Nfor j=1:nu

12、mex=ex+x(i)-rex(j,i);ey=ey+y(i)-rey(j,i);endex1(i)=ex/num;ey1(i)=ey/num;ex=0;eqx=0;ey=0;eqy=0;end%绘图figure(1);plot(x,y,k-,zx,zy,g:,xks,yks,r-.);legend(真实轨迹,观测样本,估计轨迹);figure(2);plot(ey1);legend(x方向平均误差);=基于LMS和RLS算法的自适应FIR滤波器仿真一、自适应滤波原理自适应滤波器是指利用前一时刻的结果，自动调节当前时刻的滤波器参数，以适应信号和噪声未知或随机变化的特性，得到有效的输出，主要由参

13、数可调的数字滤波器和自适应算法两部分组成，如图1.1所示自适应滤波器原理图x(n)称为输入信号，y(n)称为输出信号，d(n)称为期望信号或者训练信号，e(n)为误差僖号，其中，e(n)=d(n)-y(n)，自适应滤波器的系数(权值)根据误差信号e(n)，通过一定的自适应算法不断的进行更新，以达到使滤波器实际输出y(n)与期望响应d(n)之间的均方误差最小。1、RLS的Matlab实现：混有噪声的signalnoise.wav和噪声信号noise.wav。经过滤波之后得到filtersignal.wav既是滤波后的信号。程序如下%RLS算法，用噪声训练滤波器权系数矢量，再进行滤波%读取混合有噪

14、声的语音信号xinit,Fs,nbits = wavread(signalnoise.wav);xinit = xinit; %转置xn = xinit(1,:); %取混有噪声的语音信号的单道n = wavread(noise.wav); %读取噪声信号n = n; %转置n = n(1,:); %取噪声信号的单声道N=length(n); %噪声信号的长度M=32; %滤波器的阶数L=0.98; %遗忘因子Ns = length(xn); %有噪声的语音信号d = 0.002*ones(1,N); %期望输出x1 = 10*n+d; %RLS滤波器输入信号w=zeros(1,M); %权系

15、数x=zeros(1,M); %1*M，代表x(n)u=zeros(1,N); %临时数组u=x(n)*T(n-1)*x(n)T=200*eye(M); %T(-1)k=zeros(M,1); %增益k(n)e=zeros(1,N); %e(n|n-1)e1=e; %平方误差w1 = zeros(1,Ns-M+1);for n=M:N x=x1(n:-1:n-M+1); %x(n) u(n)=x*T*x; k=T*x/(L+u(n); %计算增益k(n) e(n)=d(n)-x*w; %计算e(n|n-1) w=w+k*e(n); %计算权系数w(n) w1(N-n+1)=w(1); T=(T

16、-x*k*T)/L; %计算T(n) e1(n)=e(n)*e(n);ends = zeros(1,Ns-M+1); %利用训练的滤波器系数进行语音信号的滤波for n=M:Ns s(n-M+1) = 0; for j=1:M s(n-M+1)= s(n-M+1)+w(j)*xn(n-j+1); endendsubplot(2,1,1); %绘制原始语音信号plot(xn);title(混有噪声的原始语音信号);subplot(2,1,2); %绘制滤波后的信号plot(s,k);title(滤波后语音信号);wavwrite(s,Fs,nbits,filter.wav);wavwrite(5*s,Fs,nbits,filtersignal.wav);下图为原始信号和除噪之后的信号对比RLS滤波前后的噪声频谱和噪声频谱比对2、LMS的Matlab实现：跟RLS的混合噪声信号相同，滤波后的噪声为signalLMS.wav%本程序只能对WAV格式的波形文件进行处理f=wavread(signalnoise.wav); a=size(f);u=0.002; %收敛步长M=64; %滤波器阶数N=a(1,1);X=z