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