MATLAB在數(shù)字信號(hào)處理中的應(yīng)用

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

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

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

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

5、設(shè)計(jì)帶通濾波器設(shè)計(jì)第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 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(相對(duì)頻率);ylab

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

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

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

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

10、為Fs=2000Hz。4.2 IIR 4.2 IIR 帶阻濾波器設(shè)計(jì)帶阻濾波器設(shè)計(jì)第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用4程序代碼1) 心電圖信號(hào)濾波代碼 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; %存有高頻干擾心電圖信號(hào)序列 k=1; close all; Figure(1) subplot

11、(2,2,1) axis(0 56 -100 50);4.2 IIR 4.2 IIR 帶阻濾波器設(shè)計(jì)帶阻濾波器設(shè)計(jì)第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用n=0:55;stem(n,x,.);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel(n);ylabel(x(n);title(心電圖信號(hào)采樣序列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 %查找行號(hào) end for l=1:4; w

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

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

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

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

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

17、的比較直接型與級(jí)聯(lián)型濾波器的比較第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 改變一對(duì)零、極點(diǎn)的位置，相對(duì)于直接型結(jié)構(gòu)，其優(yōu)點(diǎn)是調(diào)整方便，此外，級(jí)聯(lián)結(jié)構(gòu)中后面的網(wǎng)絡(luò)輸出不會(huì)流到前面，運(yùn)算誤差的積累相對(duì)于直接型也小。 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 直接型與級(jí)聯(lián)型濾波器的比較直接型與級(jí)聯(lián)型濾波器的比較第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用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(級(jí)聯(lián)型的階躍響應(yīng)); hold off;subplot(2,1,2)stem(hdir);hold on;plot(0 6,0,0);axis(0 6 -1.8 1);title(直接型的階躍響應(yīng)); hold off;第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用2) stepseq 函數(shù) function x,n=stepseq(n0,n1,n2) %Generate x(n)=u(n-n0

19、);n1=n n0;3) casfilter 函數(shù) 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在數(shù)字信號(hào)處理中的應(yīng)用4) dir2cas 函數(shù) 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在數(shù)字信號(hào)處理中的應(yīng)用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在數(shù)字信號(hào)處理中的應(yīng)用 5運(yùn)行結(jié)果與分析 b0 = 0.0667 B = 1 0 0 1 0 0 A = 1 0 0 1 0 0 由圖4.18可知，這兩種形式結(jié)構(gòu)的濾波器本質(zhì)上是完全相同的，但相比于直接型結(jié)構(gòu)，級(jí)聯(lián)形式的濾波器可以更靈活地調(diào)整零、極點(diǎn)的系數(shù)來(lái)改變一對(duì)零、極點(diǎn)的位置，調(diào)整方便，而且級(jí)聯(lián)結(jié)構(gòu)中后面的網(wǎng)絡(luò)輸出不會(huì)流到前面，運(yùn)算誤差的積累也很小，因而得到了廣泛應(yīng)用。第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 圖4.18 級(jí)聯(lián)型濾波器

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

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

24、應(yīng)濾波自適應(yīng)濾波LMSLMS算法算法第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 均方誤差是權(quán)向量的二次函數(shù)，它是一個(gè)上凹的拋物面，具有唯一的最小值，調(diào)解權(quán)向量使得均方誤差最小，相當(dāng)于沿拋物面下降尋找最小值。可以用梯度法來(lái)求該最小值，對(duì)權(quán)向量W 求導(dǎo)得到均方誤差的梯度為： 在性能曲面上最佳權(quán)矢量對(duì)應(yīng)點(diǎn)的梯度等于零，即：4.9 4.9 自適應(yīng)濾波自適應(yīng)濾波LMSLMS算法算法第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 該方程稱為正則方程，由此解出最佳權(quán)向量稱為維納解。 利用上式求解，需要精確地知道輸入信號(hào)和期望信號(hào)的先驗(yàn)統(tǒng)計(jì)知識(shí)，而且還要對(duì)矩陣求逆運(yùn)算。最陡下降法可避免求逆運(yùn)算，它通過(guò)遞推的方式尋

25、求加權(quán)矢量的最優(yōu)值，是LMS算法的理論基礎(chǔ)。首先設(shè)置一個(gè)W的初值W(0)，可以想象，沿減小的方向 調(diào)整W，可以找到最佳權(quán)矢量。因?yàn)樘荻确较蚴窃黾幼羁斓姆较?，所以?fù)梯度方向就是減少最快的方向。 最小均方算法是一種很有用很簡(jiǎn)單的估計(jì)梯度的方法，其突出特點(diǎn)是計(jì)算量小、易于實(shí)現(xiàn)，且不要求脫線計(jì)算。 LMS最核心的算法是使用平方誤差代替均方誤差，即4.9 4.9 自適應(yīng)濾波自適應(yīng)濾波LMSLMS算法算法第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用所以將上式代入最陡下降法迭代計(jì)算權(quán)矢量的公式得：則LMS算法的基本關(guān)系式為： 3仿真思路 設(shè)計(jì)一個(gè)二階加權(quán)自適應(yīng)橫向?yàn)V波器，對(duì)一個(gè)正弦信號(hào)進(jìn)行濾波。實(shí)驗(yàn)通過(guò)設(shè)置不

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

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

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

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

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

31、 頻率采樣定理頻率采樣定理第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 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);%設(shè)定頻率向量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 點(diǎn)頻率采樣幅度譜)x32n=ifft(X32k);nx32=0:31;subplot(3,2,4);stem(nx32,x32n);title(32 點(diǎn)IDFT)第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 X16k=fft(x,16); k2=0:15; subplot(3,2,5); stem(k2,abs(X16k); axis(0,8,0,200); title(16 點(diǎn)頻率采樣幅度譜) x16n=ifft(X16k); nx16=0:15; subplot(3,2,6); stem(nx16,x16n); axis(0,2

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

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

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

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

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

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

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

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

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

42、xx(500-i+j); endend第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 xd=cos(sita);%產(chǎn)生維納濾波中x 方向上觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)矩陣 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方向上最小均方誤差%產(chǎn)生維納濾波中y 方向上觀測(cè)信號(hào)的自相關(guān)矩陣 ryy=xcorr(

43、y);第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用 for i=1:100 for j=1:100 mryy(i,j)=ryy(500-i+j); end end yd=sin(sita);%產(chǎn)生維納濾波中y方向上觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)矩陣 ryd=xcorr(y,yd); for i=1:100 mryd(i)=ryd(499+i); end hopty=inv(mryy)*mryd; %由維納-霍夫方程得到的y方向上的濾波器最優(yōu)解 第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用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 方向期望信號(hào));subplot(2,4,2)plot(xnoise);title(x 方向噪聲信號(hào));subplot(2,4,3)plot(x);title(x 方向觀測(cè)信號(hào));subplot(2,4,4)n=0:500;plot(n,eminx);第4章 MATLAB在數(shù)字信號(hào)處理中的應(yīng)用title(x 方向最小均方誤差);subplot(2,4,5)plot(yd);title(y 方向期望信號(hào));subplot(2,4,6)plot(ynoise);title(y 方向噪聲信號(hào));su

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

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

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

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