维纳、卡尔曼滤波简介及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_tem

3、p; figure(1;subplot(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,ims

4、how(Blurred,title('生成的运动的模糊的图像'noise = 0.1*randn(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 =

5、 sum(noise(:.2/sum(im2double(I(:.2;wnr2 = deconvwnr(BlurredNoisy,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'

6、BW4 = edge(I,'roberts' figure(3subplot(2,4,3 4 7 8,imshow(My_Wedge,title('应用维纳滤波进行边沿提取'subplot(241,imshow(BW1,title('prewitt'subplot(242,imshow(BW2,title('canny'subplot(245,imshow(BW3,title('zerocross'subplot(246,imshow(BW4,title('roberts'%*维纳滤波应用于图像增强

7、*for i = 1 2 3 4 5K = wiener2(I,5,5;endK = K + I;figure(4;subplot(121,imshow(I,title('原始图像'subplot(122,imshow(K,title('增强后的图像'卡尔曼滤波卡尔曼滤波的一个典型实例是从一组有限的，对物体位置的，包含噪声的观察序列预测出物体的坐标位置及速度。卡尔曼滤波利用目标的动态信息,设法去掉噪声的影响，得到一个关于目标位置的好的估计。这个估计可以是对当前目标位置的估计(滤波，也可以是对于将来位置的估计(预测，也可以是对过去位置的估计(插值或平滑。状态估计是

8、卡尔曼滤波的重要组成部分。一般来说，根据观测数据对随机量进行定量推断就是估计问题，特别是对动态行为的状态估计，它能实现实时运行状态的估计和预测功能。受噪声干扰的状态量是个随机量，不可能测得精确值，但可对它进行一系列观测，并依据一组观测值，按某种统计观点对它进行估计。使估计值尽可能准确地接近真实值，这就是最优估计。真实值与估计值之差称为估计误差。卡尔曼提出的递推最优估计理论，采用状态空间描述法，在算法采用递推形式，卡尔曼滤波能处理多维和非平稳的随机过程。卡尔曼滤波的matlab实现如下所示：%卡尔曼滤波器的应用，运动轨迹估计clc;clear;T=4;%雷达扫描周期num=100;%滤波次数%产

9、生真实轨迹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+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(

10、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 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 滤波开始f

11、or 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=(eye(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 %结束一次

12、滤波endex=0;ey=0;eqx=0;eqy=0;ey1=0;ex1=N:1;ey1=N:1;%计算滤波的均值，计算滤波误差的均值for i=1:Nfor j=1:numex=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('真实轨迹','观测样本','估计轨迹

13、9;figure(2;plot(ey1;legend('x方向平均误差'=基于LMS和RLS算法的自适应FIR滤波器仿真一、自适应滤波原理自适应滤波器是指利用前一时刻的结果，自动调节当前时刻的滤波器参数，以适应信号和噪声未知或随机变化的特性，得到有效的输出，主要由参数可调的数字滤波器和自适应算法两部分组成，如图1.1所示自适应滤波器原理图x(n称为输入信号，y(n称为输出信号，d(n称为期望信号或者训练信号，e(n为误差僖号，其中，e(n=d(n-y(n，自适应滤波器的系数(权值根据误差信号e(n，通过一定的自适应算法不断的进行更新，以达到使滤波器实际输出y(n与期望响应d(n

14、之间的均方误差最小。1、RLS的Matlab实现：混有噪声的signalnoise.wav和噪声信号noise.wav。经过滤波之后得到filtersignal.wav既是滤波后的信号。程序如下%RLS算法，用噪声训练滤波器权系数矢量，再进行滤波%读取混合有噪声的语音信号xinit,Fs,nbits = wavread('signalnoise.wav'xinit = xinit' %转置xn = xinit(1,:; %取混有噪声的语音信号的单道n = wavread('noise.wav' %读取噪声信号n = n' %转置n = n(1,:

15、; %取噪声信号的单声道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; %权系数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(-1k=zeros(M,1; %增益k(ne=zeros(1,N; %e(n|n-1e1=e; %平方误差w1 = zeros(1,Ns-

16、M+1;for n=M:Nx=x1(n:-1:n-M+1; %x(nu(n=x*T*x'k=T*x'/(L+u(n; %计算增益k(ne(n=d(n-x*w' %计算e(n|n-1w=w+k'*e(n; %计算权系数w(nw1(N-n+1=w(1;T=(T-x*k*T/L; %计算T(ne1(n=e(n*e(n;ends = zeros(1,Ns-M+1; %利用训练的滤波器系数进行语音信号的滤波for n=M:Nss(n-M+1 = 0;for j=1:Ms(n-M+1= s(n-M+1+w(j*xn(n-j+1;endendsubplot(2,1,1; %绘

17、制原始语音信号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=zeros(M,1;W=zeros(M,1;out=ze