數(shù)字信號處理實驗報告 完整版

./實驗1利用DFT分析信號頻譜一、實驗?zāi)康?.加深對DFT原理的理解。2.應(yīng)用DFT分析信號的頻譜。3.深刻理解利用DFT分析信號頻譜的原理,分析實現(xiàn)過程中出現(xiàn)的現(xiàn)象及解決方法。二、實驗設(shè)備與環(huán)境計算機、MATLAB軟件環(huán)境三、實驗基礎(chǔ)理論1.DFT與DTFT的關(guān)系有限長序列的離散時間傅里葉變換在頻率區(qū)間的N個等間隔分布的點上的N個取樣值可以由下式表示：由上式可知,序列的N點DFT,實際上就是序列的DTFT在N個等間隔頻率點上樣本QUOTE。2.利用DFT求DTFT方法1：由QUOTE恢復(fù)出QUOTE的方法如下：IDFTDTFTIDFTDTFT??<??????>QUOTE由圖2.1所示流程可知：由上式可以得到：其中QUOTE為插函數(shù)方法2：實際在MATLAB計算中,上述插值運算不見得是最好的辦法。由于DFT是DTFT的取樣值,其相鄰兩個頻率樣本點的間距為2π/N,所以如果我們增加數(shù)據(jù)的長度N,使得到的DFT譜線就更加精細,其包絡(luò)就越接近DTFT的結(jié)果,這樣就可以利用DFT計算DTFT。如果沒有更多的數(shù)據(jù),可以通過補零來增加數(shù)據(jù)長度。3.利用DFT分析連續(xù)信號的頻譜采用計算機分析連續(xù)時間信號的頻譜,第一步就是把連續(xù)信號離散化,這里需要進行兩個操作：一是采樣,二是截斷。對于連續(xù)時間非周期信號QUOTE,按采樣間隔T進行采樣,階段長度M,那么：對QUOTE進行N點頻域采樣,得到因此,可以將利用DFT分析連續(xù)非周期信號頻譜的步驟歸納如下：〔1確定時域采樣間隔T,得到離散序列〔2確定截取長度M,得到M點離散序列,這里為窗函數(shù)?！?確定頻域采樣點數(shù)N,要求N≥M?！?利用FFT計算離散序列的N點DFT,得到.〔5根據(jù)上式由計算采樣點QUOTE的近似值。采用上述方法計算信號QUOTE的頻譜需要注意如下三個問題：〔1頻譜混疊。如果不滿足采樣定理的條件,頻譜會出現(xiàn)混疊誤差。對于頻譜無限寬的信號,應(yīng)考慮覆蓋大部分主要頻率分量的圍。〔2柵欄效應(yīng)和頻譜分辨率。使用DFT計算頻譜,得到的結(jié)果只是N個頻譜樣本值,樣本值之間的頻譜是未知的,像通過一個柵欄觀察頻譜,稱為"柵欄效應(yīng)"。頻譜分辨率與記錄長度成反比,要提高頻譜分辨率,就要增加記錄時間。〔3頻譜泄露。對信號截斷會把窗函數(shù)的頻譜引入信號頻譜,造成頻譜泄露。解決這個問題的主要辦法是采用旁瓣小的窗函數(shù),頻譜泄露和窗函數(shù)均會引起誤差。因此,要合理選取采樣間隔和截取長度,必要時還需考慮加適當(dāng)?shù)拇啊τ谶B續(xù)時間周期信號,我們在采用計算機進行計算時,也總是要進行截斷,序列總是有限長的,仍然可以采用上述方法近似計算。4.可能用到的MATLAB函數(shù)與代碼實驗中DFT運算可采用MATLAB中提供的函數(shù)fft來實現(xiàn)。DTFT可采用MATLAB矩陣運算的方法進行計算,如下式所示：四、實驗容1、已知,完成如下要求：〔1計算其DTFT,并畫出區(qū)間的波形。〔2計算4點DFT,并把結(jié)果顯示在〔1所畫的圖形中。〔3對x<n>補零,計算64點DFT,并顯示結(jié)果?！?根據(jù)實驗結(jié)果,分析是否可以由DFT計算DTFT,如果可以,如何實現(xiàn)。解：〔1計算其DTFT,并畫出區(qū)間的波形。>>n=0:3;>>x=[2-111];>>w=-pi:0.01*pi:pi;>>X=x*exp<-j*n'*w>;>>subplot<211>;>>plot<w,abs<X>>;>>xlabel<'\Omega/\pi'>;>>title<'Magnitude'>;>>axistight;>>subplot<212>;>>plot<w,angle<X>/pi>;>>xlabel<'\Omega/\pi'>;>>title<'Phase'>;>>axistight;〔2計算4點DFT,并把結(jié)果顯示在〔1所畫的圖形中。>>Xk=fft<x>;>>subplot<211>;>>holdon;>>stem<n,abs<Xk>,'filled'>;>>plot<w,abs<X>>;>>axistight;>>xlabel<'\Omega/\pi'>;>>title<'Magnitude'>;>>subplot<212>;>>holdon;>>plot<w,angle<X>/pi>;>>stem<n,angle<Xk>,'filled'>;>>axistight;>>xlabel<'\Omega/\pi'>;>>title<'Phase'>;運行結(jié)果如下：〔3對x<n>補零,計算64點DFT,并顯示結(jié)果。>>x=[2-111zeros<1,60>];>>n=0:63;>>Xk=fft<x>;>>subplot<211>;>>stem<n,abs<Xk>,'filled'>;>>subplot<212>;>>stem<n,angle<Xk>,'filled'>;〔4根據(jù)實驗結(jié)果,分析是否可以由DFT計算DTFT,如果可以,如何實現(xiàn)。 可以由DFT計算DTFT,序列補零后,長度越長,DFT點越多,其DFT越逼近DTFT連續(xù)波形。所以令序列補零至足夠長時其DFT序列的波形與DTFT的波形在一定的分辨率下已經(jīng)相同。2、考察序列x<n>=cos<0.48πn>+cos<0.52πn>〔10<=n<=10時,用DFT估計x<n>的頻譜；將x<n>補零加長到長度為100點序列用DFT估計x<n>的頻譜,要求畫出相應(yīng)波形。〔20<=n<=100時,用DFT估計x<n>的頻譜。并畫出波形?！?根據(jù)實驗結(jié)果,分析怎樣提高頻譜分辨率解：〔10<=n<=10時,用DFT估計x<n>的頻譜；將x<n>補零加長到長度為100點序列用DFT估計x<n>的頻譜,要求畫出相應(yīng)波形。>>n=0:10;>>x=cos<0.48*pi.*n>+cos<0.52*pi.*n>;>>Xk=fft<x>;>>subplot<211>;>>stem<n,Xk,'filled'>;>>x=[cos<0.48*pi.*n>+cos<0.52*pi.*n>zeros<1,89>];>>Xk=fft<x>;>>subplot<212>;>>stem<Xk,'filled'>;.〔20<=n<=100時,用DFT估計x<n>的頻譜。并畫出波形。程序代碼如下：>>n=0:100;>>x=cos<0.48*pi.*n>+cos<0.52*pi.*n>;>>Xk=fft<x>;>>stem<Xk,'filled'>;>>axistight;〔3根據(jù)實驗結(jié)果,分析怎樣提高頻譜分辨率可以通過如下三種方式來增加分辨率。a、增加時域信號采樣時間b、提高采樣頻率c、補零3、已知信號x<t>=0.15sin<2πf1t>+sin<2πf2t>-0.1sin<2πf3t>,其f1=1Hz,f2=2Hz,f3=3Hz。從x<t>的表達式可以看出,它包含三個頻率的正弦波,但是,從其時域波形來看,似乎是一個正弦信號,利用DFT做頻譜分析,確定適合的參數(shù),使得到的頻譜的頻率分辨率符合需要。n=0:10;x=0.15*sin<2*pi.*n>+sin<4*pi.*n>-0.1*sin<6*pi.*n>;Xk=fft<x>;stem<abs<Xk>,'filled'>;n=0:100;x=0.15*sin<0.2*pi.*n>+sin<0.4*pi.*n>-0.1*sin<0.6*pi.*n>;Xk=fft<x>;stem<abs<Xk>,'filled'>;n=0:200;x=0.15*sin<0.1*pi.*n>+sin<0.2*pi.*n>-0.1*sin<0.3*pi.*n>;Xk=fft<x>;stem<abs<Xk>,'filled'>;結(jié)果分析：上圖為x<t>信號截取過后的連續(xù)時間信號的傅里葉變換幅頻特性曲線,截取周期為100s〔即為采樣時間M,根據(jù)的特點,知1hz處的幅值為,,與圖像相符。4、利用DFT近似分析連續(xù)時間信號x<t>=e-0.1u<t>的頻譜〔幅度值。分析采用不同的采樣間隔和截取長度進行計算的結(jié)果,并最終確定合適的參數(shù)。n=0:10;x=exp<-0.1.*n>;Xk=fft<x>;stem<Xk,'filled'>;n=0:20;x=exp<-0.05.*n>;Xk=fft<x>;stem<Xk,'filled'>;n=0:40;x=exp<-0.025.*n>;Xk=fft<x>;stem<Xk,'filled'>;五、心得體會通過本次實驗,加深了對DFT原理的理解。學(xué)會了應(yīng)用DFT分析信號的頻譜。深刻理解到利用DFT分析信號頻譜的原理,能夠分析實現(xiàn)過程中出現(xiàn)的現(xiàn)象及解決方法。實驗二利用FFT計算線性卷積一、實驗?zāi)康?.掌握利用FFT計算線性卷積的原理及具體實現(xiàn)方法。2.加深理解重疊相加法和重疊保留法。3.考察利用FFT計算線性卷積各種方法的適用圍。二、實驗基礎(chǔ)理論1.線性卷積與圓周卷積設(shè)x<n>為L點序列,h<n>為M點序列,x<n>和h<n>的線性卷積為的長度為L+M-1。x<n>和h<n>的圓周卷積為圓周卷積與線性卷積相等而不產(chǎn)生交疊的必要條件為N圓周卷積定理：根據(jù)DFT性質(zhì),x<n>和h<n>的N點圓周卷積的DFT等于它們的DFT的乘積：2.快速卷積快速卷積發(fā)運用圓周卷積實現(xiàn)線性卷積,根據(jù)圓周卷積定理利用FFT算法實現(xiàn)圓周卷積。可將快速卷積運算的步驟歸納如下：<1>必須選擇；為了能使用基-2算法,要求。采用補零的辦法使得x<n>和h<n>的長度均為N。<2>計算x<n>和h<n>的N點FFT。<3>組成乘積<3>利用IFFT計算Y<k>的IDFT,得到線性卷積y<n>3.分段卷積我們考察單位取樣響應(yīng)為h<n>的線性系統(tǒng),輸入為x<n>,輸出為y<n>,則如果x<n>極長時,如果要等x<n>全部集齊時再開始進行卷積,會使輸出有較大延時；如果序列太長,需要大量存儲單元。為此,我們把x<n>分段,為別求出每段的卷積,合在一起得到最后的總輸出。這稱為分段卷積。分段卷積可以細分為重疊保留法和重疊相加法。重疊保留法：設(shè)x<n>的長度為,h<n>的長度為M。把序列x<n>分成多段N點序列,每段雨前一段重寫M-1個樣本。并在第一個輸入段前面補M-1個零。計算每一段與h<n>的圓周卷積,其結(jié)果中前M-1個不等與線性卷積,應(yīng)當(dāng)舍去,只保留后面N-M+1個正確的輸出樣本,把它們合起來得到總的輸出。利用FFT實現(xiàn)重疊保留法的步驟如下：<1>在x<n>前面填充M-1個零,擴大以后的序列為<2>將x<n>分為若干段N點子段,設(shè)L=N-M+1為每一段的有效長度,則第i段的數(shù)據(jù)為：<3>計算每一段與h<n>的N點圓周卷積,利用FFT計算圓周卷積<4>設(shè)每一段卷積結(jié)果的前M-1個樣本,連接剩下的樣本得到卷積結(jié)果y<n>。重疊相加法：設(shè)h<n>長度為M,將信號x<n>分解成長為L的子段。以表示沒斷信號,則：每一段卷積的長度為L+M-1,所以在做求和時,相鄰兩段序列由M-1個樣本重疊,即前一段的最后M-1個樣本和下一段前M-1個樣本序列重疊,這個重疊部分相加,再與不重疊的部分共同組成y<n>。利用FFT實現(xiàn)重疊保留法的步驟如下：<1>將x<n>分為若干L點子段。<2>計算每一段與h<n>的卷積,根據(jù)快速卷積法利用FFT計算卷積。<3>將各段相加,得到輸出y<n>。三、實驗容假設(shè)要計算序列和的線性卷積,完成以下實驗容：設(shè)L=M,根據(jù)線性卷積的表達式和快速卷積的原理,分別編程實現(xiàn)計算兩個序列線性卷積的方法,比較當(dāng)序列長度分別為8,16,32,64,256,512,1024時,兩種計算方法計算線性卷積所需時間。>>fori=1:7L=input<'L:'>;n=0:L;x=heaviside<n>-heaviside<n-L>;h=cos<0.2*pi.*n>;ticy=conv<x,h>;tocendL:8時間已過0.000050秒。L:16時間已過0.000037秒。L:32時間已過0.000056秒。L:64時間已過0.000065秒。L:256時間已過0.000101秒。L:512時間已過0.000188秒。L:1024時間已過0.000254秒。>>fori=1:7L=input<'L:'>;n=0:L;x=heaviside<n>-heaviside<n-L>;h=cos<0.2*pi.*n>;ticXk=fft<x,L+1>;Hk=fft<h,L+1>;Yk=Xk.*Hk;y=ifft<Yk>;tocendL:8時間已過0.000048秒。L:16時間已過0.000044秒。L:32時間已過0.000042秒。L:64時間已過0.000055秒。L:256時間已過0.000134秒。L:512時間已過0.000147秒。L:1024時間已過0.000176秒。當(dāng)L=2048且M=256時,比較計算線性卷積和快速卷積所需的時間,進一步考察當(dāng)L=4096且M=256時兩種算法所需時間。>>fori=1:2L=input<'L:'>;M=input<'M:'>;n1=0:L;n2=0:M;x=heaviside<n1>-heaviside<n1-L>;h=cos<0.2*pi.*n2>;ticXk=fft<x,L+1>;Hk=fft<h,L+1>;Yk=Xk.*Hk;y=ifft<Yk>;tocendL:2048M:256時間已過0.001139秒。L:4096M:256時間已過0.001503秒。>>fori=1:2L=input<'L:'>;M=input<'M:'>;n1=0:L;n2=0:M;x=heaviside<n1>-heaviside<n1-L>;h=cos<0.2*pi.*n2>;ticy=conv<x,h>;tocendL:2048M:256時間已過0.000263秒。L:4096M:256時間已過0.000368秒。3.>>L=input<'L:'>;M=input<'M:'>;n2=0:M;k=L/M;h=cos<0.2*pi.*n2>;Y=0;ticfori=1:kn1=<i-1>*M:i*M-1;x=heaviside<n1>-heaviside<n1-L>;Xk=fft<x,M>;Hk=fft<h,M>;Yk=Xk.*Hk;y=ifft<Yk,M>;Y=Y+y;endtocL:2048M:256時間已過0.025870秒。>>L=input<'L:'>;M=input<'M:'>;n2=0:M;k=L/M;h=cos<0.2*pi.*n2>;Y=0;ticfori=1:kn1=<i-1>*M:i*M-1;x=heaviside<n1>-heaviside<n1-L>;Xk=fft<x,M>;Hk=fft<h,M>;Yk=Xk.*Hk;y=ifft<Yk,M>;Y=Y+y;endtocL:4096M:256時間已過0.022300秒。L=4096;n1=0:L;x=heaviside<n1>-heaviside<n1-L>;n2=0:256;h=cos<0.2*pi.*n2>;ticy1=conv<x,h>;tocticXk=fft<x>;Hk=fft<h>;Yk=Xk.*Hk;y2=ifft<Yk>;toc編程實現(xiàn)利用重疊相加法計算兩個序列的線性卷積,考察L=2048且M=256時計算線性卷積的時間,與第二題結(jié)果進行比較。clearticN=512;m=0:256;h=cos<0.2*pi*m>;n=0:2048;x=heaviside<n>-heaviside<n-2048>;Lenx=length<x>;M=length<h>;M1=M-1;L=N-M1;h=fft<h,N>;K=ceil<Lenx/L>;fori=Lenx:K*L-1x<i+1>=0;endY=zeros<K,N>;YY=zeros<1,<K-1>*L+N>;fork=0:K-1xk=[x<k*L+1:k*L+L>,zeros<1,M1>];Y<k+1,:>=<ifft<fft<xk>.*h>>;YY<k*L+1:k*L+N>=YY<k*L+1:k*L+N>+Y<k+1,:>;endtoc時間已過0.028816秒。編程實現(xiàn)利用重疊保留法計算兩個序列的線性卷積,考察L=2048且M=256時計算線性卷積的時間,與第二題結(jié)果進行比較。clcclearticN=512;m=0:256;h=cos<0.2*pi*m>;n=0:2048;x=heaviside<n>-heaviside<n-2048>;Lenx=length<x>;M=length<h>;M1=M-1;L=N-M1;h=fft<h,N>;K=floor<<Lenx+M1-1>/L>+1;p=<K>*L-Lenx;x1=[zeros<1,M1>,x,zeros<1,p>];Y=zeros<K,N>;fork=0:K-1xk=fft<x1<k*L+1:k*L+N>>;Y<k+1,:>=<ifft<xk.*h>>;endZ=reshape<Y<:,M:N>',1,[]>;toc時間已過0.044424秒。四、實驗心得通過本次實驗,掌握了利用FFT計算線性卷積的原理及具體實現(xiàn)方法,加深了理解重疊相加法和重疊保留法。實驗3IIR數(shù)字濾波器設(shè)計一、實驗?zāi)康恼莆绽妹}沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的原理及具體方法。加深理解數(shù)字濾波器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化。掌握脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用圍。二、實驗容設(shè)采樣頻率為10kHz,設(shè)計數(shù)字低通濾波器,滿足如下指標(biāo) 通帶截止頻率：1kHz,通帶波動：1dB： 阻帶截止頻率：1.5kHz,阻帶衰減：15dB：要求分別采用巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和橢圓模擬原型濾波器,并分別結(jié)合脈沖響應(yīng)不變法和雙線性變換法進行設(shè)計。結(jié)合實驗結(jié)果,分別討論采用上述方法設(shè)計的數(shù)字濾波器是否都能滿足給定指標(biāo)要求,分析脈沖響應(yīng)不變法和雙線性變換法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用圍。1巴特沃斯脈沖響應(yīng)不變法：fp=1000fs=1500f=10000Wp=2*pi*fpWs=2*pi*fswp=2*f*tan<2*pi*fp/<2*f>>ws=2*f*tan<2*pi*fs/<2*f>>Rp=1As=15N1=ceil<<log10<<10^<Rp/10>-1>/<10^<As/10>-1>>>/<2*log10<Wp/Ws>>>N2=ceil<<log10<<10^<Rp/10>-1>/<10^<As/10>-1>>>/<2*log10<wp/ws>>>Omegac1=Wp/<<10^<Rp/10>-1>^<1/<2*N1>>>Omegac2=wp/<<10^<Rp/10>-1>^<1/<2*N1>>>[z,p,k]=buttap<N1>b1=k*Omegac1^N1a1=poly<p*Omegac1>[H,w]=freqs<b1,a1>subplot<221>plot<w/pi,abs<H>>gridonsubplot<223>plot<w/pi,angle<H>>gridon[b2,a2]=butter<N2,Omegac2,'S'>[H,w]=freqs<b2,a2>subplot<222>plot<w/pi,abs<H>>axistightgridonsubplot<224>plot<w/pi,angle<H>>gridon雙線性法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;N=ceil<<1/2>*<log10<<10^<As/10>-1>/<10^<Rp/10>-1>>>/<log10<ws1/wp1>>>wc=wp1/<10^<Rp/10>-1>^<1/2/N>[b1,a1]=butter<N,wc,'low','s'>;[b,a]=bilinear<b1,a1,F>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;2切比雪夫I型脈沖響應(yīng)不變法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=wp/T;ws1=ws/T;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=wp/pi[b,a]=cheby1<N,Rp,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;雙線性變換法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=wp/pi[b,a]=cheby1<N,Rp,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;3切比雪夫II型脈沖響應(yīng)不變法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=wp/T;ws1=ws/T;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=ws/pi[b,a]=cheby2<N,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;雙線性變換法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=ws/pi[b,a]=cheby2<N,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;4橢圓原型模擬濾波器脈沖響應(yīng)不變法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;wp1=wp/T;ws1=ws/T;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;k=wp1/ws1;k1=e/<A^2-1>^<1/2>;K1=ellipke<k>;K2=ellipke<<1-k1^2>^<1/2>>K3=ellipke<k1>K4=ellipke<<1-k^2>^<1/2>>;N=ceil<K1*K2/<K3*K4>>wc=wp/pi[b,a]=ellip<N,Rp,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;雙線性法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;k=wp1/ws1;k1=e/<A^2-1>^<1/2>;K1=ellipke<k>;K2=ellipke<<1-k1^2>^<1/2>>K3=ellipke<k1>K4=ellipke<<1-k^2>^<1/2>>;N=ceil<K1*K2/<K3*K4>>wc=wp/pi[b,a]=ellip<N,Rp,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'|H<e^j^\Omega>|'>;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'|H<e^j^\Omega>|<dB>'>;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'PhaseofH<e^j^Omega><\pi>'>;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'Groupdelay'>;三、結(jié)果分析：〔1對于設(shè)計要求的符合情況 根據(jù)以dB為單位的幅頻特性圖可以看出,能滿足設(shè)計要求的濾波器有：脈沖響應(yīng)不變法的切比雪夫Ⅰ型,以及雙線性變換法的全部四種濾波器。而脈沖響應(yīng)不變法的巴特沃斯、切比雪夫Ⅱ型和橢圓數(shù)字濾波器不滿足設(shè)計要求?！?脈沖響應(yīng)不變法的優(yōu)缺點及適用圍優(yōu)點：幅頻特性與模擬濾波器線性對應(yīng)。脈沖響應(yīng)不變法的一個重要特點是頻率坐標(biāo)的變換是線性的,ω=ΩΤ,ω與Ω是線性關(guān)系。因此如果模擬濾波的頻響帶限于折疊頻率以的話,通過變換后濾波器的頻響可不失真地反映原響應(yīng)與頻率的關(guān)系。穩(wěn)定性好。如果Ha<s>是穩(wěn)定的,即其極點在S左半平面,映射到H<Z>也是穩(wěn)定的。缺點：設(shè)計結(jié)果不固定。在設(shè)計時,往往不能直接根據(jù)要求的參數(shù)得到合適的濾波器,需要設(shè)計者多次修改設(shè)計參數(shù)。可能產(chǎn)生混疊失真。因頻譜周期延拓效應(yīng),脈沖響應(yīng)不變法只能用于帶限的頻響特性,如衰減特性很好的低通或帶通。而高頻衰減越大,頻響的混淆效應(yīng)越小,至于高通和帶限濾波器,由于它們在高頻部分不衰減,因此將完全混淆在低頻響應(yīng)中。所以用脈沖響應(yīng)不變法實現(xiàn)高通和帶限濾波器時,應(yīng)增加一保護濾波器,濾掉高于折疊頻率以上的頻帶,然后再用脈沖響應(yīng)不變法轉(zhuǎn)換為數(shù)字濾波器,這會增加設(shè)計的復(fù)雜性和濾波器的階數(shù)。適用圍：根據(jù)幅頻特性與模擬濾波器相對應(yīng)的特點,在要求時域脈沖響應(yīng)能模仿模擬濾波器的場合,一般使用脈沖響應(yīng)不變法?！?雙線性變換法的優(yōu)缺點及適用圍優(yōu)點：消除了脈沖響應(yīng)不變法所固有的混疊誤差?？款l率的嚴重非線性關(guān)系得到S平面與Z平面的單值一一對應(yīng)關(guān)系,整個jΩ軸單值對應(yīng)于單位圓一周,其中ω和Ω為非線性關(guān)系。設(shè)計結(jié)果具有一般性。雙線性變換法根據(jù)設(shè)計要求和計算流程計算得到的濾波器,一般能夠滿足原來的要求,得到的頻率特性較好。雙線性變換比脈沖響應(yīng)法的設(shè)計計算更直接和簡單。由于s與z之間的簡單代數(shù)關(guān)系,所以從模擬傳遞函數(shù)可直接通過代數(shù)置換得到數(shù)字濾波器的傳遞函數(shù)。缺點：頻率響應(yīng)有畸變。Ω與ω的非線性關(guān)系,導(dǎo)致數(shù)字濾波器的幅頻響應(yīng)相對于模擬濾波器的幅頻響應(yīng)有畸變,<使數(shù)字濾波器與模擬濾波器在響應(yīng)與頻率的對應(yīng)關(guān)系上發(fā)生畸變>。例如,一個模擬微分器,它的幅度與頻率是線性關(guān)系,但通過雙線性變換后,就不可能得到數(shù)字微分器。一個線性相位的模擬濾波器經(jīng)雙線性變換后,濾波器就不再有線性相位特性。適用圍：雙線性變換目前仍是使用得最普遍、最有成效的一種設(shè)計工具。這是因為大多數(shù)濾波器都具有分段常數(shù)的頻響特性,如低通、高通、帶通和帶阻等,它們在通帶要求逼近一個衰減為零的常數(shù)特性,在阻帶部分要求逼近一個衰減為∞的常數(shù)特性,這種特性的濾波器通過雙線性變換后,雖然頻率發(fā)生了非線性變化,但其幅頻特性仍保持分段常數(shù)的特性。四、實驗心得通過本次實驗,學(xué)會了掌握利用脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的原理及具體方法。加深了理解數(shù)字濾波器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化。掌握到脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的優(yōu)缺點及適用圍實驗4FIR數(shù)字濾波器設(shè)計一、實驗?zāi)康恼莆沾昂瘮?shù)法設(shè)計FIR數(shù)字濾波器的原理及具體方法。二、實驗設(shè)備與環(huán)境計算機、MATLAB軟件環(huán)境。三、實驗容1.設(shè)計一個數(shù)字低通濾波器FIR濾波器,其技術(shù)指標(biāo)如下,,分別采用矩形窗、漢寧窗、海明窗、布萊克曼窗、凱瑟窗設(shè)計該濾波器。結(jié)合實驗結(jié)果,分別討論采用上述方法設(shè)計的數(shù)字濾波器是否都能滿足給定指標(biāo)要求。〔1矩形窗程序代碼如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%錄入給定的參數(shù)tr_width=wst-wp;N=ceil<1.8*pi/tr_width>+1%根據(jù)性能指標(biāo)確定濾波器長度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%確定理想濾波器的單位脈沖響應(yīng)hd<n>w_boxcar=boxcar<N>';h=hd.*w_boxcar;%產(chǎn)生矩形窗函數(shù),并對hd<n>加窗得到h<n>subplot<221>;%畫出時域和頻域圖形,以確定設(shè)計出來的濾波器是否符合給定的指標(biāo)stem<n,hd,'filled'>;[Hr,wl]=zerophase<h>;subplot<222>;plot<wl/pi,Hr>;axistight;xlabel<'\omega/\pi'>;ylabel<'H<\omega>'>;subplot<223>;stem<n,h,'filled'>;axistight;xlabel<'n'>;ylabel<'h<n>'>;[H,w]=freqz<h,1>subplot<224>;plot<w/pi,20*log10<abs<H>/max<abs<H>>>>;xlabel<'\omega/\pi'>;ylabel<'dB'>;gridon;%set<gca,'xtick',<0:0.2:1>>;line<[0.3,0.3],[0-100],'linestyle',':'>;line<[0,1],[-0.25-0.25],'linestyle',':'>;line<[0,1],[-50-50],'linestyle',':'>;%axis<[01-20]>;運行結(jié)果如下：N=19〔2漢寧窗關(guān)鍵代碼如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%錄入給定的參數(shù)tr_width=wst-wp;N=ceil<6.2*pi/tr_width>+1%根據(jù)漢寧窗的性能指標(biāo)確定濾波器的長度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%確定理想濾波器的單位脈沖響應(yīng)hd<n>w_hanning=hanning<N>';h=hd.*w_hanning;%產(chǎn)生矩形窗函數(shù),并對hd<n>加窗得到h<n>運行結(jié)果如下：N=63〔3海明窗關(guān)鍵代碼如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%錄入給定的參數(shù)tr_width=wst-wp;N=ceil<6.6*pi/tr_width>+1%根據(jù)海明窗的性能指標(biāo)確定濾波器的長度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%確定理想濾波器的單位脈沖響應(yīng)hd<n>w_hamming=hamming<N>';h=hd.*w_hamming;%產(chǎn)生矩形窗函數(shù),并對hd<n>加窗得到h<n>實驗結(jié)果如下：N=67〔4布萊克曼窗關(guān)鍵代碼如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%錄入給定的參數(shù)tr_width=wst-wp;N=ceil<11*pi/tr_width>+1%根據(jù)布萊克曼的性能指標(biāo)確定濾波器的長度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%確定理想濾波器的單位脈沖響應(yīng)hd<n>w_blackman=blackman<N>';h=hd.*w_blackman;%產(chǎn)生矩形窗函數(shù),并對hd<n>加窗得到h<n>運行結(jié)果如下：N=111〔5凱瑟窗wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%錄入給定的參數(shù)tr_width=wst-wp;N=ceil<<As-7.95>/<2.285*tr_width>>+1%根據(jù)凱瑟窗的性能指標(biāo)確定濾波器的長度Nbeta=0.1102*<As-8.7>;n=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%確定理想濾波器的單位脈沖響應(yīng)hd<n>w_kaiser=kaiser<N,beta>';h=hd.*w_kaiser;%產(chǎn)生矩形窗函數(shù),并對hd<n>加窗得到h<n>運行結(jié)果如下：N=60上面5種設(shè)計方案中,矩形窗和漢寧窗不滿足設(shè)計要求,而海明窗、布萊克曼窗和凱瑟窗滿足設(shè)計要求。2.設(shè)計一個數(shù)字帶通FIR濾波器,其技術(shù)指標(biāo)如下下阻帶邊緣：,下阻帶邊緣：,上通帶邊緣：,上阻帶邊緣：,用布萊克曼窗設(shè)計FIR數(shù)字濾波器程序代碼如下：wp1=0.2*pi;Rp=1;wst1=0.35*pi;As=60;%錄入給定的參數(shù)wp2=0.65*pi;wst2=0.8*pi;tr_width1=wst1-wp1;tr_width2=wst2-wp2;N1=ceil<11*pi/tr_width1>+1;%根據(jù)性能指標(biāo)確定濾波器長度NN2=ceil<11*pi/tr_width2>+1;N=max<N1,N2>%帶通濾波器長度N應(yīng)為N1、N2中的較大值n=0:<N-1>;wc1=<wp1+wst1>/2;wc2=<wp2+wst2>/2;alpha=<N-1>/2;hd=<wc2/pi>*sinc<<wc2/pi>*<n-alpha>>-<wc1/pi>*sinc<<wc1/pi>*<n-alpha>>w_blackman=blackman<N>';h=hd.*w_blackman;%產(chǎn)生矩形窗函數(shù),并對hd<n>加窗得到h<n>subplot<221>;stem<n,hd,'filled'>;[Hr,wl]=zerophase<h>;subplot<222>;plot<wl/pi,Hr>;axistight;xlabel<'\omega/\pi'>;ylabel<'H<\omega>'>;subplot<223>;stem<n,h,'filled'>;axistight;xlabel<'n'>;ylabel<'h<n>'>;[H,w]=freqz<h,1>;subplot<224>;plot<w/pi,20*log10<abs<H>/max<abs<H>>>>;xlabel<'\omega/\pi'>;ylabel<'dB'>;gridon;set<gca,'xtick',<0:0.2:1>>;line<[0.35,0.35],[0-100],'linestyle',':'>;line<[0.65,0.65],[0-100],'linestyle',':'>;line<[0,1],[-1-1],'linestyle',':'>;line<[0,1],[-60-60],'linestyle',':'>;axis<[0.30.7-20]>;運行結(jié)果如下：N=753.采用頻率取樣設(shè)計法設(shè)計FIR數(shù)字低通濾波器,滿足以下指標(biāo)〔1取N=20,過渡帶沒有樣本?！?取N=40,過渡帶有一個樣本,T=0.39?！?取N=60,過渡帶有兩個樣本,T1=0.5925,T2=0.1099?！?分別討論采用上述方法設(shè)計的數(shù)字濾波器是否都能滿足給定指標(biāo)要求。解：〔1取N=20,過渡帶沒有樣本。關(guān)鍵代碼如下：N=20;alpha=<N-1>/2;l=0:N-1;wl=<2*pi/N>*l;Hrs=[1,1,1,zeros<1,15>,1,1];Hdr=[1100];wdl=[00.250.251];%用于繪制理想濾波器幅度函數(shù)曲線k1=0:floor<N-1>/2;k2=floor<<N-1>/2>+1:N-1;angH=[<-alpha*<2*pi>/N*k1>,alpha*<2*pi>/N*<N-k2>];H=Hrs.*exp<angH*j>;h=ifft<H,N>;w=[0:500]*pi/500;H=freqz<h,1,w>;[Hrwr]=zerophase<h>;運行結(jié)果如下：〔2取N=40,過渡帶有一個樣本,T=0.39。關(guān)鍵代碼如下：N=40;alpha=<N-1>/2;l=0:N-1;wl=<2*pi/N>*l;Hrs=[ones<1,5>,0.39,zeros<1,29>,0.39,ones<1,4>];Hdr=[1100];wdl=[00.250.251];%用于繪制理想濾波器幅度函數(shù)曲線k1=0:floor<N-1>/2;k2=floor<<N-1>/2>+1:N-1;angH=[<-alpha*<2*pi>/N*k1>,alpha*<2*pi>/N*<N-k2>];H=Hrs.*exp<angH*j>;h=ifft<H,