MATLAB在數字信號處理中的應用

1、第4章 MATLAB在數字信號處理中的應用 MATLAB 中的信號處理工具箱提供了很多數字信號處理中需要用到的函數及解決方法。在數字信號處理中所學習到的問題，諸如濾波器的設計、自適應濾波，維納濾波、卡爾曼濾波等理論都可以通過MATLAB仿真得到實現(xiàn)和驗證。本章提供了12個MATLAB在數字信號處理課程學習中的應用實例，可作為學習數字信號處理課程的參考。第4章 MATLAB在數字信號處理中的應用知知 識識 架架 構構 1實驗目的 (1)掌握脈沖響應不變法設計IIR數字濾波器的具體設計方法。 (2)熟悉脈沖響應不變法設計低通濾波器的仿真。 2實驗原理 脈沖響應不變法是從濾波器的脈沖響應出發(fā)，使數字

2、濾波器的單位脈沖響應序列h(n)模仿模擬濾波器的沖擊響應ha(t)，使h(n)正好等于ha(t)的采樣值，即：4.1 IIR 4.1 IIR 帶通濾波器設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用 T為采樣周期。若以Ha(s)及H(z)分別表示ha(t)的拉氏變換及h(n)的Z變換，即：4.1 IIR 4.1 IIR 帶通濾波器設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用 根據采樣序列Z變換與模擬信號拉氏變換的關系得: 上式表明，采用脈沖響應不變法將模擬濾波器變換為數字濾波器時，它所完成的S平面到Z平面的變換，正是以前討論的拉氏變換到Z變換的標準變換關系。 脈沖

3、響應不變法特別適用于用部分分式表達的傳遞函數，模擬濾波器的傳遞函數若只有單階極點，且分母的階數高于分子階數，則可表達為部分分式形式：4.1 IIR 4.1 IIR 帶通濾波器設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用 其拉氏反變換為： 對ha(t)采樣就得到數字濾波器的單位脈沖響應序列： 再對h(n)取Z變換，得到數字濾波器的傳遞函數： 第二個求和為等比級數之和：4.1 IIR 4.1 IIR 帶通濾波器設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用 要收斂的話，必有： 比較部分分式形式的Ha(s)和上式H(z)可以看到，把S平面上的極點變換到Z平面上對應的極點

4、，而Ha(s)與H(z)中部分分式所對應的系數不變。如果模擬濾波器是穩(wěn)定的，則所有極點都在S左半平面那么變換后H(z)的極點也都在單位圓內，因此數字濾波器保持穩(wěn)定。 所以有： 3仿真思路 在MATLAB中，可以用下列函數輔助設計IIR數字濾波器。 (1)利用cheb1ord可以確定低通原型和切比雪夫濾波器的階數和截止頻率。 (2)num，den=cheby1(N，Wn)，num，den=cheby2(N，Wn) (切比雪夫1 型和2 型)可以進行濾波器的設計。 (3)利用impinvar可以完成脈沖響應不變法的模擬濾波器到數字濾波器的轉換。 4程序代碼4.1 IIR 4.1 IIR 帶通濾波器

5、設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用 clear wp=6*pi*103;ws=9*pi*103;ap=1,as=15; Fs=30*103; wp1=wp/Fs;ws1=ws/Fs; N,WC=cheb1ord(wp,ws,ap,as,s); b,a=cheby1(N,ap,WC,s) ； bz,az=impinvar(b,a,Fs)； w0=wp1,ws1； Hx=freqz(bz,az,w0)； H,W=freqz(bz,az)； dbHx=-20*log10(abs(Hx)/max(abs(H)； plot(W,abs(H); xlabel(相對頻率);ylab

6、el(幅頻); grid4.1 IIR 4.1 IIR 帶通濾波器設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用 5運行結果與分析 結論：由圖4.2可知，切比雪夫型濾波器的振幅特性在通帶內是等波紋的，在阻帶內是單調的。4.1 IIR 4.1 IIR 帶通濾波器設計帶通濾波器設計第4章 MATLAB在數字信號處理中的應用6思考題(1) 帶通濾波器的幅度有什么特點？圖4.2中濾波器的下降沿怎么樣才能更陡峭？(2) 用本節(jié)所學的知識設計切比雪夫型濾波器。圖 4.2： 1實驗目的 (1)熟悉用雙線性變換法設計IIR數字濾波器的原理與方法。 (2)熟悉帶阻數字濾波器設計方法。 (3)通過觀

7、察對實際心電圖信號的濾波作用，獲得數字濾波的感性知識。 2實驗原理 1)用雙線性變換法設計IIR數字濾波器數字低通技術指標為：4.2 IIR 4.2 IIR 帶阻濾波器設計帶阻濾波器設計第4章 MATLAB在數字信號處理中的應用 模擬低通的技術指標為：查巴特沃斯歸一化低通濾波參數表可得：第4章 MATLAB在數字信號處理中的應用設計巴特沃斯低通濾波器。階數N計算如下。數字濾波器系統(tǒng)函數H(z)為：4.2 IIR 4.2 IIR 帶阻濾波器設計帶阻濾波器設計 A=0.09036 B1=1.2686, C1=0.7051 B2=1.0106, C2=0.3583 B3=0.9044, C3=0.2

8、155 可見H(z)是由3個二階濾波器H1(z)、H2(z)、H3(z)級聯(lián)組成的，如圖4.3所示。4.2 IIR 4.2 IIR 帶阻濾波器設計帶阻濾波器設計第4章 MATLAB在數字信號處理中的應用 圖 4.3 濾波器H(z)的組成圖2) 帶阻IIR數字濾波器實現(xiàn)過程 (1)按一定規(guī)格將數字濾波器的技術指標轉為模擬低通濾波器的技術指標。 (2)根據轉換后的技術指標使用濾波器階數函數，確定最小階數N 和截止頻率Wc。 (3)利用最小階數N產生模擬低通濾波器原型。 (4)利用截止頻率Wc把模擬低通原型轉化為模擬帶阻濾波器。 (5)利用沖激響應不變法或雙線性不變法把模擬濾波器轉換成數字濾波器。

9、3仿真思路 (1)用雙線性變換法設計一個低通IIR數字濾波器。4.2 IIR 4.2 IIR 帶阻濾波器設計帶阻濾波器設計第4章 MATLAB在數字信號處理中的應用 設計參數為：在通帶內頻率低于0.2rad時，最大衰減小于1dB，在阻帶內0.3，頻率區(qū)間上，最小衰減大于15dB。 (2)所設計的濾波器對實際心電圖信號采樣序列進行仿真濾波處理，觀察濾波前后的圖形。 (3)設計一個帶阻IIR數字濾波器，其具體要求是：通帶的截止頻率wp1=650Hz，wp2=850Hz；阻帶的截止頻率ws1=700Hz，ws2=800Hz，帶通內的最大衰減rp=0.1dB，阻帶內的最小衰減為rs=50dB，采樣頻率

10、為Fs=2000Hz。4.2 IIR 4.2 IIR 帶阻濾波器設計帶阻濾波器設計第4章 MATLAB在數字信號處理中的應用4程序代碼1) 心電圖信號濾波代碼 x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,. -4,-4,-6,-6,-2, 6,12, 8, 0,-16,. -38,-60,-84,-90,-66,-32,-4,-2,-4, 8,. 12,12,10,6,6,6,4,0,0,0,. 0,0,-2,-4,0,0,0,-2,-2,0,. 0,-2,-2,-2,-2,0; %存有高頻干擾心電圖信號序列 k=1; close all; Figure(1) subplot

11、(2,2,1) axis(0 56 -100 50);4.2 IIR 4.2 IIR 帶阻濾波器設計帶阻濾波器設計第4章 MATLAB在數字信號處理中的應用n=0:55;stem(n,x,.);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel(n);ylabel(x(n);title(心電圖信號采樣序列x(n);B=0.09036 2*0.09036 0.09036;A=1.2686 -0.7051;A1=1.0106 -0.3583;A2=0.9044 -0.2155;while(k0.7,break,end %查找行號 end for l=1:4; w

12、pz2=0.4 0.3625 0.325 0.3; %低頻的4 個高通截止數字頻率 wsz2=0.25; rp2=3;rs2=60;第4章 MATLAB在數字信號處理中的應用 N2,wc2=buttord(wpz2(l),wsz2,rp2,rs2); %高通數字濾波器 Bz2,Az2=butter(N2,wc2,high); y2=filter(Bz2,Az2,x); %濾波后的單頻信號 z2=max(abs(y2); %濾波后信號在離散域所取的最大值 subplot(3,3,m+6);plot(n,y2);grid;title( 高通后的波形);xlabel(n); ylabel(y) if

13、 z20.7,break,end; %查找列號 end TNr=TNr+tm(k,(5-l)*10(3-m); %將3位電話號碼表示成一個3位數 end disp(接受斷檢測到的號碼為：) disp(TNr) 5運行結果與分析 從圖4.17中可知，3個數字的低頻部分頻譜相同， 高頻部分不同，跟假設匹配。通過對濾出頻率分析，查找表，可以得出輸出的數字。第4章 MATLAB在數字信號處理中的應用 圖4.17數字1、2、3的時域和濾波后的波形 6思考題 雙音多頻濾波器常用來解決什么問題？第4章 MATLAB在數字信號處理中的應用 1實驗目的 (1)知道直接型濾波器轉換為級聯(lián)型濾波器的原理。 (2)了

14、解直接型濾波器轉換為級聯(lián)型濾波器后，級聯(lián)型網絡結構的優(yōu)點。 (3)掌握直接型濾波器轉換為級聯(lián)型濾波器的實現(xiàn)方式及級聯(lián)型濾波器函數輸出的編程原理。 2實驗原理 為了對連續(xù)的或離散的信號進行濾波處理，就必須構造出合適的實際結構。對于同樣的系統(tǒng)函數H(z)或H(s)往往有不同的實現(xiàn)方案。常用的有直接形式和級聯(lián)形式以及并聯(lián)形式。下面主要討論直接形式4.8 4.8 直接型與級聯(lián)型濾波器的比較直接型與級聯(lián)型濾波器的比較第4章 MATLAB在數字信號處理中的應用 到級聯(lián)形式的轉換以及二者的比較。 直接形式用延遲元件和乘法器以及加法器以給定的形式直接實現(xiàn)差分方程，設M=N=4，則差分方程為：4.8 4.8 直

15、接型與級聯(lián)型濾波器的比較直接型與級聯(lián)型濾波器的比較第4章 MATLAB在數字信號處理中的應用 在這種形式中將系統(tǒng)函數H(z)寫成實系數二階子系統(tǒng)的乘積形式。首先把分子和分母多項式的根解出，然后把每一對共軛復根或任意兩個實根組合在一起，得到二階子系統(tǒng)。在以下原理中假設N 為偶 數，于是可把上式轉化為下面的形式：4.8 4.8 直接型與級聯(lián)型濾波器的比較直接型與級聯(lián)型濾波器的比較第4章 MATLAB在數字信號處理中的應用上式中的參量滿足：二階子系統(tǒng)為： 在工程實際中，一般把上式所示的結構稱為雙二階環(huán)節(jié)(Biquad),它的輸入是第(k 1) 個雙二階環(huán)節(jié)的輸出，同時第k個雙二階環(huán)節(jié)的輸出為第(k

16、+1)個雙二階環(huán)節(jié)的輸入。 3仿真思路 本程序由1個主函數和3個子函數組成，主函數完成由下面的差分方程描述的直接型濾波器到級聯(lián)型濾波器的轉換。 15y(n)+10y(n1)+ 3y(n2)5y(n3)2y(n4) =x(n)4x(n1)+10 x(n2)25x(n3) +16x(n4) 子函數dir2cas實現(xiàn)將直接形式轉化為級聯(lián)形式，子函數stepseq(n0,n1,n2)，實現(xiàn)階躍函數u(n-n0)，n1nn2。子函數casfilter是編制實現(xiàn)級聯(lián)形式濾波器的函數。 級聯(lián)結構中每一個二階網絡決定一對零點和一對極點，可以通過靈活地調整對應零、極點的系數來4.8 4.8 直接型與級聯(lián)型濾波器

17、的比較直接型與級聯(lián)型濾波器的比較第4章 MATLAB在數字信號處理中的應用 改變一對零、極點的位置，相對于直接型結構，其優(yōu)點是調整方便，此外，級聯(lián)結構中后面的網絡輸出不會流到前面，運算誤差的積累相對于直接型也小。 4程序代碼 1)主程序 b=1 -4 10 -25 16; a=15 10 3 -5 -2; b0,B,A=dir2cas(b,a) stepin=stepseq(0,0,5); format long; hcas=casfilter(b0,B,A,stepin)4.8 4.8 直接型與級聯(lián)型濾波器的比較直接型與級聯(lián)型濾波器的比較第4章 MATLAB在數字信號處理中的應用hdir=f

18、ilter(b,a,stepin)format;figure(1)subplot(2,1,1)stem(hcas); hold onplot(0 6,0 0);axis(0 6 -1.8 1);title(級聯(lián)型的階躍響應); hold off;subplot(2,1,2)stem(hdir);hold on;plot(0 6,0,0);axis(0 6 -1.8 1);title(直接型的階躍響應); hold off;第4章 MATLAB在數字信號處理中的應用2) stepseq 函數 function x,n=stepseq(n0,n1,n2) %Generate x(n)=u(n-n0

19、);n1=n n0;3) casfilter 函數 function y=casfilter(b0,B,A,x) K,L=size(B); N=length(x); w=zeros(K+1,N); w(1,:)=x; for i=1:K w(i+1,:)=filter(B(i,:),A(i,:),w(i,:); end y=b0*w(K+1,:);第4章 MATLAB在數字信號處理中的應用4) dir2cas 函數 function b0,B,A=dir2cas(b,a) b0=b(1); b=b0/b0; a0=a(1); a=a/a0; b0=b0/a0; M=length(b); N=l

20、ength(a); if NM b=b zeros(1,N-M); else if NM a=a zeros(1,N-M); else NM=0; end K=floor(N/2);第4章 MATLAB在數字信號處理中的應用B=zeros(K,3);A=zeros(K,3);if K*2=M b=b 0; a=a 0;endbroots=cplxpair(roots(b);aroots=cplxpair(roots(a);For i=1:2:2*K; Brow=broots(i:1:i+1,:); Brow=real(poly(Brow); B(fix(i+1)/2),:)=Brow; Aro

21、w=aroots(i:1:i+1,:); Arow=real(poly(Arow); A(fix(i+1)/2),:)=Brow;End第4章 MATLAB在數字信號處理中的應用 5運行結果與分析 b0 = 0.0667 B = 1 0 0 1 0 0 A = 1 0 0 1 0 0 由圖4.18可知，這兩種形式結構的濾波器本質上是完全相同的，但相比于直接型結構，級聯(lián)形式的濾波器可以更靈活地調整零、極點的系數來改變一對零、極點的位置，調整方便，而且級聯(lián)結構中后面的網絡輸出不會流到前面，運算誤差的積累也很小，因而得到了廣泛應用。第4章 MATLAB在數字信號處理中的應用 圖4.18 級聯(lián)型濾波器

22、和直接型濾波器的階躍響應 6思考題 (1)級聯(lián)型濾波器和直接型濾波器各有什么優(yōu)缺點？ (2)級聯(lián)型濾波器是如何實現(xiàn)的？第4章 MATLAB在數字信號處理中的應用 1實驗目的 (1)進一步了解自適應濾波原理。 (2)學習LMS自適應算法及其MATLAB仿真。 2實驗原理 自適應濾波器由參數可調的數字濾波器和自適應算法兩部分組成。輸入信號x(n)通過參數可調數字濾波器后產生輸出信號y(n)，將其與參數信號d(n)進行比較，形成誤差信號e(n)。e(n)通過某種自適應算法對濾波器參數進行調整，最終使e(n)的均方值最小。 最小均方誤差LMS準則的目的在于使濾波器輸出與期望信號誤差的平方的統(tǒng)計平均值最

23、小。LMS自適應橫向濾波器如圖4.19所示。4.9 4.9 自適應濾波自適應濾波LMSLMS算法算法第4章 MATLAB在數字信號處理中的應用 圖4.19 最小均方誤差準則的自適應橫向濾波器原理圖該自適應濾波器的輸入矢量為：4.9 4.9 自適應濾波自適應濾波LMSLMS算法算法第4章 MATLAB在數字信號處理中的應用加權矢量為為：濾波器的輸出為：y(n)相對于濾波器期望輸出d(n)的誤差為： 根據最小均方誤差準則，最佳的濾波器參量應使得性能函數均方誤差為最小。 假定輸入信號x(n)和期望響應d(n)是聯(lián)合平穩(wěn)過程，那么在時刻n的均方誤差是加權矢量的二次函數，其表示式為：4.9 4.9 自適

24、應濾波自適應濾波LMSLMS算法算法第4章 MATLAB在數字信號處理中的應用 均方誤差是權向量的二次函數，它是一個上凹的拋物面，具有唯一的最小值，調解權向量使得均方誤差最小，相當于沿拋物面下降尋找最小值?？梢杂锰荻确▉砬笤撟钚≈?，對權向量W 求導得到均方誤差的梯度為： 在性能曲面上最佳權矢量對應點的梯度等于零，即：4.9 4.9 自適應濾波自適應濾波LMSLMS算法算法第4章 MATLAB在數字信號處理中的應用 該方程稱為正則方程，由此解出最佳權向量稱為維納解。 利用上式求解，需要精確地知道輸入信號和期望信號的先驗統(tǒng)計知識，而且還要對矩陣求逆運算。最陡下降法可避免求逆運算，它通過遞推的方式尋

25、求加權矢量的最優(yōu)值，是LMS算法的理論基礎。首先設置一個W的初值W(0)，可以想象，沿減小的方向 調整W，可以找到最佳權矢量。因為梯度方向是增加最快的方向，所以負梯度方向就是減少最快的方向。 最小均方算法是一種很有用很簡單的估計梯度的方法，其突出特點是計算量小、易于實現(xiàn)，且不要求脫線計算。 LMS最核心的算法是使用平方誤差代替均方誤差，即4.9 4.9 自適應濾波自適應濾波LMSLMS算法算法第4章 MATLAB在數字信號處理中的應用所以將上式代入最陡下降法迭代計算權矢量的公式得：則LMS算法的基本關系式為： 3仿真思路 設計一個二階加權自適應橫向濾波器，對一個正弦信號進行濾波。實驗通過設置不

26、同的收斂因子，由MATLAB程序圖形觀察濾波效果。討論的重要性。 4程序代碼 clear all ; fs =10000; t=0：1/fs ：1; sn =sin(2 * pi * t); %產生初始信號 n= randn (size (t); %產生高斯噪聲 xn = sn + n; %信號加噪聲 w = 0 0.5; %設置權初值 u=0.00026; %設置收斂因子 for i =1:length(t)-1第4章 MATLAB在數字信號處理中的應用 for i =1:length(t)-1 y (i+1)=n(i:i+1)*w; %噪聲通過濾波器輸出為y e (i+1)=xn (i+1

27、)-y(i+1); w = w+2*u*e(i+1)*n(i:i+1); %權的變化公式 end subplot (3,1,1) plot (t, xn ) %輸出信號加噪聲圖形 title (帶噪聲原始信號) grid； subplot (3,1,2) plot (t , sn) title(原始正弦信號) subplot(3,1,3) plot (t , e) title(濾波結果)第4章 MATLAB在數字信號處理中的應用 5 5運行結果與分析運行結果與分析 收斂因子決定收斂速度及穩(wěn)定性。圖4.20設置的收斂因子為0.00026，圖4.21設置的收斂因子為0.0026，圖4.22設置的收

28、斂因子為0.026。由圖可見，收斂因子的選擇對濾波器的性能有很大影響。圖4.20，濾波效果不錯，但是收斂較慢；圖4.21，收斂速度加快，但濾波效果不如圖4.20；圖4.22 收斂速度很快，但濾波效果太差。可見，選擇合適的收斂因子對于濾波性能有很大影響。收斂因子的選擇，其實就是在濾波性能及速度之間折中。 6思考題 (1)簡述LMS 算法的原理。 (2)根據上述結果，收斂因子該如何選擇？第4章 MATLAB在數字信號處理中的應用 圖4.20 =0.00026 圖4.21 =0.0026第4章 MATLAB在數字信號處理中的應用 圖4.22 =0.026 時正弦加噪信號的濾波第4章 MATLAB在數

29、字信號處理中的應用 1實驗目的 (1)了解頻率采樣定理是數字信號處理中的重要理論。 (2)掌握頻率域采樣會引起時域周期化的概念。 (3)掌握頻率域采樣定理及其對頻率采樣點數選擇的 指導作用。 2實驗原理 1) 頻率域采樣定理的要點 對信號x(n)的頻譜在0,2上等間隔采樣N點， 得到:4.10 4.10 頻率采樣定理頻率采樣定理第4章 MATLAB在數字信號處理中的應用 則N點的離散傅里葉反變換得到的序列就是原序 列x(n)以N為周期進行周期延拓后的主值序列，公式 為：4.10 4.10 頻率采樣定理頻率采樣定理第4章 MATLAB在數字信號處理中的應用 由上式可知，頻率域采樣點數N必須大于等

30、于時域離散信號的長度M(即NM)，才能使時域不產生混疊，則N點的離散傅里葉反變換得到的序列就是原序列項x(n)。 對比上面敘述的頻域采樣定理，得到一個有用的結論：頻域采樣，時域信號周期延拓。 2) 頻域采樣定理的驗證 給定信號如下： 編寫程序分別對頻譜函數在區(qū)間0，2上等間隔采樣32點和16點，得到：4.10 4.10 頻率采樣定理頻率采樣定理第4章 MATLAB在數字信號處理中的應用 3仿真思路 直接調用MATLAB函數fft計算再分別對進行32點和16點IFFT，得到：分別畫出的幅度譜，并繪圖顯示的波形，進行對比和分析。 就得到頻譜函數在0，2上的32點和16點頻率采樣。4.10 4.10

31、 頻率采樣定理頻率采樣定理第4章 MATLAB在數字信號處理中的應用 4程序代碼 M=26;N=32;n=0:26; n1=0:M/2; n2=M/2+1:M; xa=n1+1; xb=27-n2; x=xa,xb; subplot(3,2,2); stem(n,x); title(三角波序列x(n) ) w=linspace(0,pi,1000);%設定頻率向量X=x*exp(-j*n*w);subplot(3,2,1);plot(w/pi,abs(X);title(序列x(n)的連續(xù)幅度譜)X32k=fft(x,32);k1=0:31;subplot(3,2,3);stem(k1,abs(

32、X32k);axis(0,15,0,200);title(32 點頻率采樣幅度譜)x32n=ifft(X32k);nx32=0:31;subplot(3,2,4);stem(nx32,x32n);title(32 點IDFT)第4章 MATLAB在數字信號處理中的應用 X16k=fft(x,16); k2=0:15; subplot(3,2,5); stem(k2,abs(X16k); axis(0,8,0,200); title(16 點頻率采樣幅度譜) x16n=ifft(X16k); nx16=0:15; subplot(3,2,6); stem(nx16,x16n); axis(0,2

33、0,0,28); title(16 點IDFT)5運行結果與分析 由圖4.23可知，在一定的條件下，可以由頻域離第4章 MATLAB在數字信號處理中的應用 散采樣恢復原來的信號，這條件就是：如果序列x(n)的長度為M，則只有當頻域采樣點數NM時，才可由頻域采樣X(k)恢復原序列x(n)，否則將產生時域混疊現(xiàn)象。 6思考題 (1) 頻率采樣定理和時域采樣定理有什么不同？ (2) 頻率采樣定理在什么情況下能恢復原來的信號？4.10 4.10 頻率采樣定理頻率采樣定理第4章 MATLAB在數字信號處理中的應用第4章 MATLAB在數字信號處理中的應用圖 4.23 對不同信號的采樣 1實驗目的 (1)

34、了解維納濾波器的原理。 (2)了解維納濾波器的MATLAB實現(xiàn)方法。 2實驗原理 濾波技術是信號分析、處理技術的重要分支，無論是信號的獲取、傳輸，還是信號的處理和交換都離不開濾波技術，它對信號安全可靠和有效靈活地傳遞是至關重要的。信號分析檢測與處理的一個十分重要的內容就是從噪聲中提取信號，實現(xiàn)這種功能的有效手段之一是設計具有最佳線性過濾特性的濾波器，當伴有噪聲的信號通過這種濾波器時，它可以將信號盡4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用 可能精確地重現(xiàn)或對信號做出盡可能精確的估計，而對所伴隨噪聲進行最大限度的抑制。維納濾波器就是這種濾波器的典型代表

35、之一。 1) 維納濾波概述 維納(Wiener)是用來解決從噪聲中提取信號的一種濾波方法。這種線性濾波問題可以看成是一種估計問題或一種線性估計問題。 如果一個線性系統(tǒng)的單位樣本響應為h(n)，當輸入一個隨機信號x(n)，且4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用其中x(n)表示信號，v(n)表示噪聲，則輸出y(n)為： 希望x(n)通過線性系統(tǒng)h(n)后得到的y(n)盡量接近于s(n)，因此稱y(n)為s(n)的估計值，即： 則維納濾波器的輸入-輸出關系可用圖4.24表示。4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處

36、理中的應用 圖4.24 維納濾波器的輸入-輸出關系 用h(n)進行過濾問題實際上是一種統(tǒng)計估計問題。一般地，從當前的和過去的觀察值x(n)，x(n 1)， x(n 2)，估計當前的信號值稱為過濾或濾波；從過去的觀察值，估計當前的或者將來的信號值稱為外推或預測；從過去的觀察值，估計過去的 信號值稱為平滑或內插。因此維納濾波器又常常被稱為最佳線性過濾與預測或線性最優(yōu)估計。這里所謂的最佳與最優(yōu)是以最小均方誤差為準則的。 用e(n)表示信號的真實值和估計值之間的誤差，顯然e(n)可能是正值，也可能是負值，并且它是一個隨機變量。因此用它的均方誤差來表達誤差是合理的，均方誤差最小即它的平方的統(tǒng)計期望最小。

37、采用最小均方誤差準則作為最佳過濾準則的原因還在于它的理論分析比較簡單，不要求對概率的描述。 2) 維納-霍夫方程的求解4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用 為了按最小均方誤差準則來確定維納濾波器的沖激響應h(n)，令(n) 對h(j)的導數等于零，即可得維納-霍夫方程：4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用其中： 維納-霍夫方程右端的求和范圍沒有具體標明，實際上有3種情況。 有限沖激響應(FIR)維納濾波器，i從0到N1取得有限個整數值。 非因果無限沖激響應(非因果IIR)維納濾波器，i從到+取所

38、有整數值。 因果無限沖激響應(因果IIR)維納濾波器，i從0到+取正整數值。 上述3種情況下維納-霍夫方程的解法不同，這里只描述FIR維納濾波器的求解。 設濾波器沖激響應序列的長度為N，沖激響應矢量為：4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用濾波器輸入數據矢量為：則濾波器的輸出為：這樣，上述所示的維納濾波器的標準方程可寫成： 其中，P是s(n)與x(n)的互相關函數，它是一個N維列矢量；R是x(n)的自相關函數，是N階方陣。 利用求逆矩陣的方法直接求解上式，得：4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用o

39、pt表示“最佳”，這就是FIR維納濾波器的沖激響應。 3仿真思路 假設一個點目標在 x、y平面上繞單位圓做圓周運動，由于外界干擾，其運動軌跡發(fā)生了偏移。其中，x方向的干擾為均值為0、方差為0.05的高斯噪聲；y方向的干擾為均值為0、方差為0.06的高斯噪聲。 (1)產生滿足要求的x方向和y方向隨機噪聲樣本500 個。 (2)明確期望信號和觀測信號。 (3)試設計一個FIR維納濾波器，確定最佳傳遞函數并用該濾波器處理觀測信號，得到其最佳估計。 (4)分別繪制出x方向和y 方向的期望信號、噪聲信號、觀測信號、濾波后信號、最小均方誤差信號的曲線圖。 (5)在同一幅圖中繪制出期望信號、觀測信號和濾波后

40、點目標的運動軌跡。 4程序代碼4.11 4.11 維納濾波算法維納濾波算法第4章 MATLAB在數字信號處理中的應用 clear; clf; sita=0:pi/249.5:2*pi; xnoise=sqrt(0.05)*randn(1,500);%產生x軸方向噪聲 ynoise=sqrt(0.06)*randn(1,500);%產生y軸方向噪聲 x=cos(sita)+xnoise;%產生x軸方向觀測信號 y=sin(sita)+ynoise;%產生y軸方向觀測信號 %產生維納濾波中x方向上觀測信號的自相關矩陣 rxx=xcorr(x); for i=1:100 for j=1:100 mr

41、xx(i,j)=rxx(500-i+j); endend第4章 MATLAB在數字信號處理中的應用 clear; clf; sita=0:pi/249.5:2*pi; xnoise=sqrt(0.05)*randn(1,500);%產生x軸方向噪聲 ynoise=sqrt(0.06)*randn(1,500);%產生y軸方向噪聲 x=cos(sita)+xnoise;%產生x軸方向觀測信號 y=sin(sita)+ynoise;%產生y軸方向觀測信號 %產生維納濾波中x方向上觀測信號的自相關矩陣 rxx=xcorr(x); for i=1:100 for j=1:100 mrxx(i,j)=r

42、xx(500-i+j); endend第4章 MATLAB在數字信號處理中的應用 xd=cos(sita);%產生維納濾波中x 方向上觀測信號與期望信號的互相關矩陣 rxd=xcorr(x,xd); for i=1:100 mrxd(i)=rxd(499+i); end hoptx=inv(mrxx)*mrxd; %由維納-霍夫方程得到的x 、 方向上的濾波器最優(yōu)解 fx=conv(x,hoptx); %濾波后x方向上的輸出 nx=sum(abs(xd).2); eminx=nx-mrxd*hoptx; %x方向上最小均方誤差%產生維納濾波中y 方向上觀測信號的自相關矩陣 ryy=xcorr(

43、y);第4章 MATLAB在數字信號處理中的應用 for i=1:100 for j=1:100 mryy(i,j)=ryy(500-i+j); end end yd=sin(sita);%產生維納濾波中y方向上觀測信號與期望信號的互相關矩陣 ryd=xcorr(y,yd); for i=1:100 mryd(i)=ryd(499+i); end hopty=inv(mryy)*mryd; %由維納-霍夫方程得到的y方向上的濾波器最優(yōu)解 第4章 MATLAB在數字信號處理中的應用fy=conv(y,hopty); %濾波后y方向上的輸出ny=sum(abs(yd).2);eminy=ny-mr

44、yd*hopty; %y方向上最小均方誤差subplot(2,4,1)plot(xd);title(x 方向期望信號);subplot(2,4,2)plot(xnoise);title(x 方向噪聲信號);subplot(2,4,3)plot(x);title(x 方向觀測信號);subplot(2,4,4)n=0:500;plot(n,eminx);第4章 MATLAB在數字信號處理中的應用title(x 方向最小均方誤差);subplot(2,4,5)plot(yd);title(y 方向期望信號);subplot(2,4,6)plot(ynoise);title(y 方向噪聲信號);su

45、bplot(2,4,7)plot(y);title(y 方向觀測信號);subplot(2,4,8)plot(n,eminy);title(y 方向最小均方誤差);figure;plot(xd,yd,k);第4章 MATLAB在數字信號處理中的應用 hold on; plot(x,y,b:); hold on; plot(fx,fy,g-); title(最終結果); 5運行結果與分析 x軸和y軸信號、噪聲及最小均方誤差如圖4.25 所示。濾波后得到的信號與原始信號和噪聲信號的對比如圖4.26所示，濾波后的結果與期望信號還是很接近的，整體上達到了最優(yōu)濾波的效果。 6思考題 (1)維納濾波原理其

46、實是一個最優(yōu)估計的問題。簡述上述程序實現(xiàn)維納濾波的過程。第4章 MATLAB在數字信號處理中的應用第4章 MATLAB在數字信號處理中的應用 圖4.25 信號、噪聲及最小均方誤差濾 圖 4.26 最終結果 (2)上例中，估計信號與原始信號之間還有一種事實上 的誤差存在，為什么？如何使這種誤差最小化？第4章 MATLAB在數字信號處理中的應用 1實驗目的 (1)了解卡爾曼濾波器的原理。 (2)學習卡爾曼濾波器的仿真實現(xiàn)方法。 2實驗原理 卡爾曼濾波器是一個“ 最優(yōu)化自回歸數據處理算法。對于解決大部分的問題，它是最優(yōu)、效率最高甚至是最有用的。其廣泛應用已經超過30年，包括機器人導航、控制、傳感器數

47、據融合甚至在軍事方面的雷達系統(tǒng)以及導彈追蹤等。近年來更被應用于計算機圖像處理，例如，面部識別、圖像分割、圖像邊緣檢測等方面。4.12 4.12 卡爾曼濾波算法卡爾曼濾波算法第4章 MATLAB在數字信號處理中的應用 卡爾曼濾波原理卡爾曼濾波原理 首先要引入一個離散控制過程的系統(tǒng)，該系統(tǒng)可用一個線性隨機微分方程來描述： X(k)=AX(k1)+BU(k)+W(k) 再加上系統(tǒng)的測量值： Z(k)=H X(k)+V(k) 上兩式子中，X(k)是k時刻的系統(tǒng)狀態(tài)，U(k)是k 時刻對系統(tǒng)的控制量。A和B是系統(tǒng)參數，對于多模型系統(tǒng)，它們?yōu)榫仃?。Z(k)是k時刻的測量值， H是測量系統(tǒng)的參數，對于多測量

48、系統(tǒng)，H為矩陣。W(k)和V(k)分別表示過程和測量的噪聲。4.12 4.12 卡爾曼濾波算法卡爾曼濾波算法第4章 MATLAB在數字信號處理中的應用 它們被假設成高斯白噪聲，其協(xié)方差分別是Q、R這 里假設它們不隨系統(tǒng)狀態(tài)變化而變化。 由于滿足上面的條件(線性隨機微分系統(tǒng)，過程和測量都是高斯白噪聲)，卡爾曼濾波器是最優(yōu)的信息處理器。下面來估算系統(tǒng)的最優(yōu)化輸出。 首先利用系統(tǒng)的過程模型預測下一個狀態(tài)的系統(tǒng)。假設現(xiàn)在的系統(tǒng)狀態(tài)是k，根據系統(tǒng)的模型，可以基于系統(tǒng)的上一狀態(tài)而預測出現(xiàn)在狀態(tài)： X(k|k-1)=AX(k-1|k-1)+BU(k) 上式中，X(k|k-1)是利用上一個狀態(tài)預測的結果，X(k-1|k-1)是上一個狀態(tài)最優(yōu)的結果，U(k)為現(xiàn)在狀4.12 4.12 卡爾曼濾波算法卡爾曼濾波算法第4章 MATLAB在數字信號處理中的應用 態(tài)的控制量，如果沒有控制量，它可以為0。 到現(xiàn)在為止，系統(tǒng)結果已經更新了，可是對應于X(k|k-1)的協(xié)方差還沒更新。用P表示協(xié)方差： P(k|k1)=A P(k1|k1) A+Q 上式中，P(k|k1)是X(k|k1)對應的協(xié)方差，P(k1|k1)是X(k1|k1)對應的協(xié)方差，A表示A的轉置矩陣，Q是系統(tǒng)過程