數(shù)字信號處理與MATLAB實現(xiàn)

1、數(shù)字信號處理與MATLAB實現(xiàn)1 n1=ns:nf；x1=zeros(1,n0-ns),1,zeros (1,nf-n0)； %單位抽樣序列的產(chǎn)生2. subplot(2,2,4) 畫2行2列的第4個圖3. stem(n,x) %輸出離散序列，(plot連續(xù))4. 編寫子程序可調(diào)用4.1 單位抽樣序列生成函數(shù)impseq.mx,m=impseq(n0,ns,nf); %序列的起點為ns，終點為nf，在n=n0點處生成一個單位脈沖n=-5:5;x1=3*impseq(2,-5,5)-impseq(-4,-5,5)x1 = 0 -1 0 0 0 0 0 3 0 0 0n=-5:5;x1=3*imp

2、seq(2,-4,5)-impseq(-4,-5,4) %起點到終點長度要一致x1 = 0 -1 0 0 0 0 3 0 0 04.2 單位階躍序列生成函數(shù)stepseq.mx,n=stepseq(no,ns,nf) %序列的起點為ns，終點為nf，在n=n0點處生成一個單位階躍4.3 兩個信號相加的生成函數(shù)sigadd.my,n=sigadd(x1,n1,x2,n2)4.4 兩個信號相乘的生成函數(shù)sigmult.my,n=sigmult(x1,n1,x2,n2)4.5 序列移位y(n)=x(n-n0)的生成函數(shù)sigshift.my,n=sigshift(x,m,n0)4.6 序列翻褶y(n

3、)=x(-n)的生成函數(shù)sigfold.my,n=sigfold(x,n)4.7 evenodd.m函數(shù)可以將任一給定的序列x(n)分解為xe(n)和xo(n)兩部分xe,xo,m=evenodd(x,n)4.8 序列從負(fù)值開始的卷積conv_m, conv默認(rèn)從0開始function y,ny=conv_m(x,nx,h,nh)有x(n):nx1nnx2,h(n):nh1nnh2, 卷積結(jié)果序列為y(n):nx1+nh1nnx2+nh2例. 設(shè),求程序：x1=1,2,3;n1=-1:1;x2=2,4,3,5;n2=-2:1;y,n=conv_m(x1,n1,x2,n2)結(jié)果：y = 2 8

4、17 23 19 15n =-3 -2 -1 0 1 2因此計算，得Z反變換b=1;a=poly(0.9 0.9 -0.9);r,p,k=residuez(b,a)結(jié)果：r = 0.2500 0.5000 0.2500p = 0.9000 0.9000 -0.9000k =因此得到 相應(yīng)的5. H,w=freqzB,A,M計算出M個頻率點上的頻率響應(yīng)，存放于H向量中，M個頻率存放在w向量中，freqz函數(shù)自動將這M個頻率點均勻設(shè)置在頻率范圍0，之間。若缺省w和M時，函數(shù)自動選取512個頻率點計算。不帶輸出向量的freqz函數(shù)將自動繪制幅頻和相頻曲線。也可H,w=freqz(B,A);plot(

5、w/pi,abs(H); 繪出幅頻特性6. zplane(z,p)繪制出列向量z中的零點（以符號o表示）和列向量p中的極點（以符號x表示）以及參考單位圓。7. 傅里葉變換信號時域和頻域變換連續(xù)對應(yīng)著非周期，離散對應(yīng)著周期。一個域的離散必然導(dǎo)致另一個域的周期延拓8. fft和ifft：一維快速傅里葉變換和逆傅里葉變換X=fft(x,N)采用FFT算法計算序列向量X的N點DFT，缺省N時，fft函數(shù)自動按X的長度計算DFT。當(dāng)N為2的整數(shù)次冪時，fft按基2算法計算，否則用混合基算法。ifft的調(diào)用格式類似。9. fft2和ifft2: 二維快速正傅里葉變換和逆傅里葉變換（1）Y=fft2(x)數(shù)

6、據(jù)二維傅里葉變換參數(shù)X是向量fft2(x)相當(dāng)于fft(fft(x),即先對X的列做一維傅里葉變換，然后再對變換結(jié)果的行作一維傅里葉變換；若X是向量，則此傅里葉變換即變成一維傅里葉變換fft; 若X是矩陣，則是計算該矩陣的二位傅里葉變換。（2）Y=fft2(X, M, N)通過對X進(jìn)行補零或截斷，使得X成為（M*N）的矩陣。函數(shù)iff2的參數(shù)應(yīng)用與函數(shù)fft2完全相同。10. czt: 線形調(diào)頻Z變換Y=czt(x,m,w,a)此函數(shù)計算由z=a*w.(-(0:m-1)定義的z平面螺旋線上各點的z變換，a規(guī)定了起點，w規(guī)定了相鄰點的比例，m規(guī)定了變換的長度，后三個變量默認(rèn)值為a=1,w=exp

7、(j*2*pi/m)及m=length(x)因此y=czt(x)就等于y=fft(x).11. dct和idct：離散余弦正變換和離散余弦逆變換Y=dct(x,N)完成如下變換，N的默認(rèn)值為length(x). k=0, 1, ，N-112. fftshiftY=fftshift(x)用來重新排列X=fft(x)的輸出，當(dāng)X為向量時，把X的左右兩半進(jìn)行交換，從而將零頻分量移至頻譜的中心；如果X為二維傅里葉變換的結(jié)果，它同時將X的左右和上下部分進(jìn)行交換。13. fftfiltY=fftfilt(b,x)采用重疊相加法FFT對信號向量x快速濾波，得到輸出序列向量y，向量b為FIR濾波的單位脈沖相應(yīng)

8、，h(n)=b(n+1),n=0,1,length(b)-1.Y=fftfilt(b,x,N)自動選取FFT長度NF=2nextpow2(N), 輸入數(shù)據(jù)x分段長度M=NF-length(b)+1, 其中nextpow2(N)函數(shù)求的一個整數(shù)，滿足2(nextpow2(N)-1)<=N<=2nextpow2(N)N缺省時，fftfilt自動選擇合適的FFT長度NF和對x的分段長度M。p = nextpow2(A)， p 滿足 2p >= abs(A)L=pow2(nextpow2(M+N-1) M, N分別為序列延拓周期14. 用DFT進(jìn)行譜分析時，必須將序列階段為長度為N的

9、有限長序列。造成頻譜泄露和譜間干擾。泄露使頻譜變得模糊，分辨率降低；旁瓣引起不同分量間的干擾。截斷效應(yīng)無法完全消除，可以加寬窗和緩慢截斷。矩形窗比海明窗的頻率分辨率高（泄漏小），但譜間干擾大，因此海明窗是以犧牲分辨率來換取譜間干擾的降低。柵欄效應(yīng)是只能在離散點的地方看到真實的像，其余頻譜被遮擋。為減少柵欄效應(yīng)，可以在時域數(shù)據(jù)末端增加一些零點，使周期內(nèi)點數(shù)增加，但不改變原有數(shù)據(jù)，即增加頻域抽樣點數(shù)N，頻域抽樣為，這樣必然使譜線更密，這樣原來看不到的譜分量就可能看到了。15. x=one(1,5),zero(1,N-5); %單位階躍信號y = filter(b,a,x) %直接型輸出信號16.

10、N=25;h=impz(b,a,N) %N次采樣，b,a為零極點多項式系數(shù)，a不算1如： %直接型b=1，-3，11，27，18；a=16，12，2，-4，-1 %級聯(lián)型b0=4;B=1,1,0;1,-1.4142136,1;A=1.-0.5.0;1.0.9.0.81 %并聯(lián)型C=0;B=-14,-12;24,26;A=1,-2,3;1,-1,1delta=impseq(0,0,N);x=one(1,5),zero(1,N-5); h=filter(b,a,delta); %直接型單位脈沖響應(yīng)y = filter(b,a,x); %直接型輸出信號h=casfilter(b0,B,A,delta

11、); %級聯(lián)型單位脈沖響應(yīng)h=casfilter(b0,B,A,x); %級聯(lián)型輸出響應(yīng)h=parfiltr(C,B,A,delta); %并聯(lián)型單位脈沖響應(yīng)h=parfiltr(C,B,A,x); %并聯(lián)型輸出響應(yīng)b0,B,A=dir2cas(b,a); %直接型轉(zhuǎn)換成級聯(lián)型b,a=cas2dir(b0,B,A); %級聯(lián)型轉(zhuǎn)換成直接型C,B,A=dir2par(b,a); %直接型轉(zhuǎn)換成并聯(lián)型b,a=par2dir(C,B,A); %并聯(lián)型轉(zhuǎn)換成直接型17. buttord.m用來確定數(shù)字低通或模擬低通濾波器的階次，其調(diào)用格式分別是(1)N,Wn=buttord(Wp,Ws,Rp,Rs)(

12、2) N,Wn=buttord(Wp,Ws,Rp,Rs,s)格式(1)對應(yīng)數(shù)字濾波器, 式中Wp,Ws分別是通帶和阻帶的截止頻率，實際上它們是歸一化頻率，其值在01之間，1對應(yīng)抽樣頻率的一半。對低通和高通濾波器，Wp,Ws都是標(biāo)量，對帶通和帶阻濾波器，Wp,Ws都是12的向量。Rp,Rs分別是通帶和阻帶的衰減，單位為dB。N是求出的相應(yīng)低通濾波器的階次，Wn是求出的3dB頻率，它和Wp稍有不同。格式（2）對應(yīng)模擬濾波器，式中各個變量的含義和格式（1）相同，但Wp,Ws及Wn的單位為rad/s, 因此，他們實際上是頻率。18. buttap.m用來設(shè)計模擬低通原型濾波器G(p), 其調(diào)用格式是z

13、,p,k=buttap(N)N是與設(shè)計的低通原型濾波器的階次，z,p和k分別是設(shè)計出的G(p)的極點、零點及增益。19. lp2lp.m, lp2hp.m, lp2bp.m, lp2bs.m以上4個文件的功能是將模擬低通原型濾波器G(p)分別轉(zhuǎn)換為實際的低通、高通、帶通及帶阻濾波器。其調(diào)用格式分別為（1）B,A=lp2lp(b,a,Wo) 或 B,A=lp2hp(b,a,Wo)（2）B,A=lp2bp(b,a,Wo,Bw) 或 B,A=lp2bs(b,a,Wo,Bw)式中b,a分別是模擬低通原型濾波器G(p)的分子、分母多項式的系數(shù)向量，B,A分別是轉(zhuǎn)換后的H(s)的分子、分母多項式的系數(shù)向量

14、；在格式（1）中，Wo是低通或高通濾波器的截止頻率；在格式（2）中Wo是帶通或帶阻濾波器的中心頻率，Bw是其帶寬。20. bilinear.m實現(xiàn)雙線性變換，即由模擬濾波器H(s)得到數(shù)字濾波器H(z)。其調(diào)用格式是Bz,Az=bilinear(B,A,Fs)式中B,A分別是H(s)的分子、分母多項式的系數(shù)向量，Bz,Az分別是H(z)的分子、分母多項式的系數(shù)向量，F(xiàn)s是抽樣頻率。21. butter.m用來直接設(shè)計巴特沃茲數(shù)字濾波器，實際上它把buttord.m, buttap.m, lp2lp.m及bilinear.m等文件都包含進(jìn)去，從而使設(shè)計過程更簡捷。其調(diào)用格式是(1) B,A=bu

15、tter(N,Wn)(2) B,A=butter(N,Wn,high)(3) B,A=butter(N,Wn,stop)(4) B,A=butter(N,Wn,s)格式(1)(3)用來設(shè)計數(shù)字濾波器，B、A分別是H(z)的分子、分母多項式的系數(shù)向量，Wn是通帶截止頻率，范圍在01之間，1對應(yīng)抽樣頻率的一半。若Wn是標(biāo)量，則格式(1)用來設(shè)計低通數(shù)字濾波器，若Wn是12的向量，則格式(1)用來設(shè)計數(shù)字帶通濾波器；格式（2）用來設(shè)計數(shù)字高通濾波器；格式（3）用來設(shè)計數(shù)字帶阻濾波器，顯然，這時的Wn是12的向量：格式（4）用來設(shè)計模擬濾波器。22. cheb1ord.m求切比雪夫1型濾波器的階次。2

16、3. cheb1ap.m用來設(shè)計原型切比雪夫1型模擬濾波器。24. cheby1.m 直接設(shè)計切比雪夫1型濾波器。以上3個文件的調(diào)用格式和對應(yīng)的巴特沃茲濾波器的文件類似。25. impinvar.m用沖激響應(yīng)不變法實現(xiàn)到及s到z的轉(zhuǎn)換。26. maxflat.m設(shè)計廣義巴特沃茲低通濾波器。27. zp2tf.m把零極點加增益的形式，改寫成零極點多項式 b,a = zp2tf(z,p,k)28. 低通數(shù)字濾波器設(shè)計例子lpexample.m模擬巴特沃茲濾波器低通濾波器設(shè)計實例monibashi.m其中，巴特沃茲模擬濾波器的設(shè)計子程序設(shè)計和非歸一化巴特沃茲模擬低通濾波器原型子程序afd_butt.

17、m計算系統(tǒng)函數(shù)的幅度響應(yīng)和相位響應(yīng)子程序freqs_m.m直接形式轉(zhuǎn)換成級聯(lián)形式子程序sdir2cas.m模擬切比雪夫1型低通濾波器設(shè)計實例moniqieshi.m（2型.m文件中改為2）其中，切比雪夫1型模擬濾波器的設(shè)計子程序設(shè)計和非歸一化切比雪夫1型模擬低通濾波器原型子程序afd_chb1.m計算系統(tǒng)函數(shù)的幅度響應(yīng)和相位響應(yīng)子程序freqs_m.m直接形式轉(zhuǎn)換成級聯(lián)形式子程序sdir2cas.m29. 用脈沖響應(yīng)不變法設(shè)計低通數(shù)字濾波器, 模擬濾波器采用切比雪夫1型濾波器原型模擬濾波器采用切比雪夫1型濾波器原型程序shuziqieshi.m脈沖響應(yīng)不變法子程序imp_invr數(shù)字濾波器響應(yīng)

18、子程序freqz_m直接型轉(zhuǎn)換成并聯(lián)型子程序dir2par.m比較兩個含同樣標(biāo)量元素,但(可能)有不同下標(biāo)的復(fù)數(shù)對及其相應(yīng)留數(shù)向量子程序cplxcomp.m用脈沖響應(yīng)不變法設(shè)計數(shù)字濾波器，模擬濾波器采用巴特沃茲濾波器原型shuzibashi.m30. 模擬低通濾波器變換為數(shù)字高通濾波器Alp2Dhp.m模擬低通濾波器變換為數(shù)字帶通濾波器Alp2Dbp.m模擬低通濾波器變換為數(shù)字帶阻濾波器Alp2Dbs.m31. 數(shù)字低通到數(shù)字濾波器的設(shè)計為了從低通數(shù)字濾波原型得到新數(shù)字濾波器的有理函數(shù)，必須實現(xiàn)有理代換，通常采用zmapping.m函數(shù)實現(xiàn)。32. 用zmapping函數(shù)實現(xiàn)高通濾波器實例zm

19、appinghp.m33. 一個含有高頻噪聲的信號保存在數(shù)據(jù)文件noisbloc中，試用脈沖響應(yīng)不變法設(shè)計低通濾波器對其進(jìn)行除噪。lpchuzao.m34. 無限長單位沖擊響應(yīng)（IIR）數(shù)字濾波器的優(yōu)點是可以利用模擬濾波器設(shè)計的結(jié)果，而模擬濾波器的設(shè)計有大量圖表可查，方便簡單。但是它也有明顯的缺點，就是相位的非線性，在圖像處理以及數(shù)據(jù)傳輸?shù)纫笮诺谰哂芯€性相位特性的場合，IIR濾波器就不太適用了。有限單位沖擊響應(yīng)（FIR）數(shù)字濾波器則可以做成具有嚴(yán)格的線性相位，同時又可以具有任意的幅度特性。此外，F(xiàn)IR濾波器的單位抽樣相應(yīng)是有限長的，因而FIR濾波器一定是穩(wěn)定的。再有，只要經(jīng)過一定的延時，任何

20、非因果有限長序列都能變成因果的有限長序列，因而總能用因果系統(tǒng)來實現(xiàn)。最后，F(xiàn)IR濾波器由于單位沖擊響應(yīng)是有限長的，因而可以用快速傅里葉變換（FFT）算法來過濾信號，從而可大大提高運算效率。但是，要取得很好的衰減特性，F(xiàn)IR濾波器H(z)的階次比IIR濾波器的要高。IIR濾波器設(shè)計中的各種變換法對FIR濾波器設(shè)計是不適用的，這是因為那里是利用有理分式的系統(tǒng)函數(shù)，而FIR濾波器的系統(tǒng)函數(shù)只是的多項式。因為最感興趣的是具有線性相位的FIR濾波器，對非線性相位的FIR濾波器，一般可以用IIR濾波器來代替。所以只討論線性相位濾波器。FIR數(shù)字濾波器的文件35. fir1.m本文件采用窗函數(shù)法設(shè)計FIR數(shù)

21、字濾波器，其調(diào)用格式是b=fir1(N,Wn)b=fir1(N,Wn,high)b=fir1(N,Wn,stop)式中N為濾波器的階次，因此濾波器的長度為N+1;Wn是通帶截止頻率，其值在01之間，1對應(yīng)抽樣頻率的一半；b是設(shè)計好的濾波器系數(shù)h(n).對于第一種格式，若Wn是一標(biāo)量，則可用來設(shè)計低通濾波器；若Wn是12的向量，則用來設(shè)計帶通濾波器；若Wn是1L的向量，則可用來設(shè)計L帶濾波器，此時，格式將變?yōu)閎=fir1(N,Wn,DC-1)或b=fir1(N,Wn,DC-0)其中，前者保證第一個帶為通帶，后者保證打一個帶為阻帶。第二種格式用來設(shè)計高通濾波器，第三種格式用來設(shè)計帶阻濾波器。值得注

22、意的是，在上述所有格式中，若不指定窗函數(shù)的類型，則fir1自動選擇漢明窗。36. fir2.m本文件采用窗函數(shù)法設(shè)計具有任意幅頻特性的FIR濾波器。其調(diào)用格式是b=fir1(N,F,M)其中F是頻率向量，其值在01之間，M是與F相對應(yīng)的所希望的幅頻響應(yīng)。不指定窗函數(shù)的類型時，將自動選擇漢明窗。37. remez.m本文件用來設(shè)計采用切比雪夫最佳一致逼近FIR數(shù)字濾波器。同時，還可以用來設(shè)計希爾伯特變換器和差分器。其調(diào)用格式是b=remez(N,F,A)b=remez(N,F,A,W)b=remez(N,F,A,W, hilbert)b=remez(N,F,A,differentiator)其中，N是給定的濾波器的階次；b是設(shè)計的濾波器的系數(shù)，其長度為N+1;F是頻率向量，其值在01之間；A是對應(yīng)F的各頻段上的理想幅頻響應(yīng)；W是各頻段上的加權(quán)向量。值得注意的是，若b的長度為偶數(shù)，涉及高通和帶阻濾波器時有可能出現(xiàn)錯誤，因此最好保證b的長度為奇數(shù)，即N應(yīng)為偶數(shù)。38. remexord.m本文件采用切比雪夫一致逼近設(shè)計FIR數(shù)字