數(shù)字信號處理matlab仿真

數(shù)字信號處理作業(yè)設(shè)計報告一、目的1.增進對Matlab的認識，加深對數(shù)字信號處理理論方面的理解。2.掌握數(shù)字信號處理中IIR和FIR濾波器的設(shè)計。3.了解和掌握用Matlab實現(xiàn)IIR和FIR濾波器的設(shè)計方法、課程，為以后的設(shè)計打下良好根底。二、設(shè)計內(nèi)容1.IIR〔無限脈沖響應(yīng)〕模擬濾波器設(shè)計〔1〕設(shè)計題目：橢圓型模擬帶通IIR濾波器技術(shù)指標：通帶下截止頻率fpl=2kHz，上截止頻率fph=5kHz，通帶內(nèi)最大衰減ap=1dB；阻帶下截止頻率fsl=1.5kHz上截止頻率fsh=5.5kHz,阻帶最小衰減as=40dB.設(shè)計原理：①確定模擬帶通濾波器的技術(shù)指標，并對邊界頻率做歸一化處理；②確定歸一化低通技術(shù)要求；③設(shè)計歸一化低通G(p)；④將低通G(p)轉(zhuǎn)換成帶通H(s)。Matlab原程序如下：clearall;fp=[2000,5000];ap=1;fs=[1500,5500];as=40;wp=2*pi*fp;ws=2*pi*fs;%歸一化的截止頻率[N,wn]=ellipord(wp,ws,ap,as,'s');%求橢圓形濾波器的最小階數(shù)和歸一化截止頻率[B,A]=ellip(N,ap,as,wn,'s');%求傳遞函數(shù)的分子分母系數(shù)[H,w]=freqs(B,A);%頻率響應(yīng)函數(shù)f=0:8000;w=2*pi*f;H=freqs(B,A,w);%求系統(tǒng)在指定頻率點w上的頻響Hplot(f,20*log10(abs(H)));%繪圖顯示axis([07000-800])仿真波形圖如下：〔2〕設(shè)計題目：巴特沃斯低通模擬濾波器技術(shù)指標：通帶截止頻率fp=5kHz，通帶內(nèi)最大衰減ap=2dB；阻帶截止頻率fs=12kHz，阻帶最小衰減as=30dB。設(shè)計原理：①確定模擬帶通濾波器的技術(shù)指標，并對邊界頻率做歸一化處理；②確定歸一化低通技術(shù)要求并求出歸一化低通原型系統(tǒng)函數(shù)Ga(p);③將Ga(p)去歸一化。Matlab原程序如下：clearall;wp=2*pi*5000;ws=2*pi*12000;Rp=2;As=30;%設(shè)置濾波器參數(shù)[N,wc]=buttord(wp,ws,Rp,As,'s');%計算濾波器階數(shù)N和3dB截止頻率wc[B,A]=butter(N,wc,'s');%計算濾波器系統(tǒng)函數(shù)分子分母多項式系數(shù)k=0:511;fk=0:14000/512:14000;wk=2*pi*fk;Hk=freqs(B,A,wk);subplot(2,2,1);plot(fk/1000,20*log10(abs(Hk)));gridonxlabel('頻率(kHz)');ylabel('幅度(dB)')axis([0,14,-40,5])仿真波形圖如下：2.IIR〔無限脈沖響應(yīng)〕數(shù)字濾波器設(shè)計〔1〕設(shè)計題目：IIR數(shù)字帶通濾波器設(shè)計技術(shù)指標：采用巴特沃斯濾波器,一個數(shù)字系統(tǒng)的抽樣頻率Fs=2000Hz，設(shè)計一個為此系統(tǒng)使用的帶通數(shù)字濾波器Hdbp(z)。要求通帶范圍為300~400Hz，在帶邊頻率處的衰減不大于3dB；在200Hz以下和500Hz以上衰減不小于18dB。設(shè)計原理：對于數(shù)字帶通濾波器的設(shè)計，通常所用方法為雙線性變換?？山柚谀M濾波器的頻率變換設(shè)計一個所需類型的過渡模擬濾波器，再通過雙線性變換將其轉(zhuǎn)換成所需類型的數(shù)字濾波器。Matlab原程序如下：clearall;fp=[300400];fs=[200500];rp=3;rs=18;Fs=2000;wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;%%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp./2);was=2*Fs*tan(ws./2);[n,wn]=buttord(wap,was,rp,rs,'s');%Note:'s'![z,p,k]=buttap(n);[bp,ap]=zp2tf(z,p,k);%bw=wap(2)-wap(1);w0=sqrt(wap(1)*wap(2));[bs,as]=lp2bp(bp,ap,w0,bw);%[h1,w1]=freqs(bp,ap);figure(1);plot(w1,abs(h1));grid;ylabel('lowpassG(p)')%w2=[0:Fs/2-1]*2*pi;h2=freqs(bs,as,w2);%Note:z=(2/Ts)(z-1)/(z+1);[bz1,az1]=bilinear(bs,as,Fs);[h3,w3]=freqz(bz1,az1,1000,Fs);figure(2);plot(w2/2/pi,20*log10(abs(h2)),w3,20*log10(abs(h3)));grid;ylabel('BandpassAFandDF')xlabel('Hz');仿真波形圖如下：圖①模擬低通原型圖圖②模擬和數(shù)字帶通濾波器的幅頻響應(yīng)〔2〕設(shè)計題目：基于Matlab函數(shù)直接設(shè)計切比雪夫Ⅰ型數(shù)字低通IIR濾波器技術(shù)指標：通帶截止頻率fp為100HZ，阻帶截止頻率fs為200HZ，通帶最大衰減ap為3dB，阻帶最小衰減as為30dB,采樣頻率Fs為1000HZ.設(shè)計原理：Matlab工具箱中提供了直接設(shè)計數(shù)字濾波器的函數(shù)使設(shè)計更加簡單。切比雪夫Ⅰ型濾波器在通帶是等波紋的，而在阻帶是單調(diào)下降的，利用cheby1函數(shù)MATLAB語法可以直接設(shè)計；而切比雪夫Ⅰ型濾波器得階次可以通過MATLAB中的cheb1ord函數(shù)確定。Matlab原程序如下：clearall;fp=100; ap=3;%通帶衰減fs=200;as=30;%阻帶衰減Fs=1000;[N,Wn]=cheb1ord(fp/(Fs/2),fs/(Fs/2),ap,as);%計算ChebyshevⅠ型濾波器的最小階數(shù)和截止頻率[b,a]=cheby1(N,ap,Wn);freqz(b,a,512,1000);%采用1000HZ的采樣頻率繪出顯示濾波器的幅頻和相頻響應(yīng)仿真波形圖如下：FIR〔有限脈沖響應(yīng)〕數(shù)字濾波器設(shè)計〔1〕設(shè)計題目：基于矩形窗和漢明窗的FIR濾波器設(shè)計技術(shù)指標：分別用矩形窗和漢明窗，設(shè)計一低通濾波器，令N=10，頻率響應(yīng)截止頻率Hd(ejω)在0≤ω≤ππ≤ω≤π之間為0。設(shè)計原理：要設(shè)計FIR濾波器，h(n)必然是有限長的，所以要用有限長的h(n)來逼近無限長的hd(n)，最有效的方法是截斷hd(n)，即用一個有限長度的窗函數(shù)ω(n)來截取hd(n)，并將截短后的hd(n)移位，得由于ω(n)的長度為N，所以h(n)是因果的。因此窗函數(shù)序列ω(n)的形狀及長度的選擇就很關(guān)鍵。Matlab原程序如下：clearall;N=10;b1=fir1(N,0.25,boxcar(N+1));%用矩形窗作為沖激響應(yīng)的窗函數(shù)b2=fir1(N,0.25,hamming(N+1));%用漢明窗作為沖激響應(yīng)的窗函數(shù)%M=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);%分別求兩個濾波器的頻率響應(yīng)t=0:10subplot(221)stem(t,b2,'.');holdon;plot(t,zeros(1,11));grid;f=0:0.5/M:0.5-0.5/M;M1=M/4;fork=1:M1hd(k)=1;hd(k+M1)=0;hd(k+2*M1)=0;hd(k+3*M1)=0;endsubplot(222)plot(f,abs(h1),'b-',f,abs(h2),'g-',f,hd,'-');grid;仿真波形圖如下：〔2〕設(shè)計題目：海明窗函數(shù)法FIR數(shù)字濾波器技術(shù)指標：ππ，阻帶最小衰減為As＝50dB。設(shè)計原理：π/N,因此具有較小的階數(shù)。在MATLAB中利用ideal_lp函數(shù)，計算理想低通濾波器的單位脈沖響應(yīng)hd〔n〕。functionhd=ideal_lp(Wc,N);alpha=(N-1)/2;n=0:1:N-1;m=n-alpha+eps;hd=sin(Wc*m)./(pi*m);利用freqz_m函數(shù)計算絕對和相對的幅度響應(yīng)。function[db,mag,pha,w]=freqz_m(b,a);[H,w]=freqz(b,a,1000,’whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);以上兩個函數(shù)應(yīng)先分別保存在兩個新建m文件中存放在Matlab軟件工作路徑下。Matlab原程序如下：clearall;Wp=0.2*pi;Ws=0.4*pi;%濾波器的邊界頻率tr_width=Ws-Wp;%過渡帶寬N=ceil(6.6*pi/tr_width)+1%根據(jù)過度帶寬等于海明窗函數(shù)主瓣寬度求得濾波器所用海明窗的最小長度n=0:1:N-1;wc=(Ws+Wp)/2;%截止頻率在通帶和阻帶邊界頻率的中點hd=ideal_lp(wc,N);%調(diào)用理想低通濾波器脈沖響應(yīng)函數(shù)w_ham=(hamming(N))';%采用海明窗h=hd.*w_ham;%在時間域乘積對應(yīng)于頻率域的卷積，截取得到實際單位脈沖響應(yīng)[db,mag,pha,w]=freqz_m(h,[1]);%調(diào)用freqz_m函數(shù)計算實際濾波器絕對和相對的幅度響應(yīng)delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1)));%實際通帶最大衰減As=-round(max(db(Ws/delta_w+1:1:501)));%實際最小阻帶衰減subplot(221);stem(n,hd);title('理想單位脈沖響應(yīng)hd(n)')subplot(222);stem(n,w_ham);title('海明窗)w(n)')subplot(223);stem(n,h);title('實際單脈沖響應(yīng)h(n)')subplot(224);plot(w/pi,db);title('幅度響應(yīng)(db)')axis([0,1,-100,10]);運行結(jié)果：N=34仿真波形圖如下：4、利用FFT進行頻譜分析設(shè)一序列中含有三種頻率成分，f1=2Hz，f2=2.05Hz，f3=1.9Hz，采樣頻率為fs=10Hz，序列：分別取N1=64,N2=128點有效數(shù)據(jù)作頻譜特性分析，分別在四個圖形窗口繪出x(n)，X(k)64點DFT，X(k)補零到128點DFT，X(k)128點DFT。比擬得出在哪種情況下可以清楚地分辨出信號的頻譜成分設(shè)計要求：因為f1=2Hz,f2=2.05Hz,f3=1.9Hz采樣頻率為fs=10Hz,所以構(gòu)造正弦序列x(n)=sin(2nπf1/fs)+sin(2nπf2/fs)+sin(2nπf3/fs).設(shè)計原理：信號采樣的頻譜分辨率為F=0.05Hz,所以采樣記錄最小點數(shù)為N=2fh/F=2f2/F=82點。利用Matlab離散傅里葉變換函數(shù)對序列進行頻譜變換，通過Matlab頻譜波形圖分析信號頻譜成分Matlab原程序如下：clearall;N1=64;fs=10;n=0:1:N1-1;f1=2;f2=2.05;f3=1.9;xn=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+sin(2*pi*f3*n/fs);%滿足條件的正弦序列x(n)Xk=fft(xn);AXk=abs(Xk(1:N1/2));%對x(n)進行64點DFT運算figure(1)subplot(2,1,1)plot(n,xn);%繪出64點采樣正弦序列x(n)figure(2)k=(0:N1/2-1)*fs/N1;subplot(2,1,2)plot(k,AXk);%繪出X(k)64點DFT頻譜幅度%X(k)補零到128點DFTN2=128;xn=[xnzeros(1,N2-N1)];%對64點采樣序列x(n)進行補零到128點Xk=fft(xn);AXk=abs(Xk(1:N2/2));%對補零到128點的序列進行DFT運算m=0:1:N2-1;figure(3)k=(0:N2/2-1)*fs/N2;subplot(2,1,2)plot(k,AXk)%繪出X(k)補零到128點DFT頻譜幅度%X(k)128點DFTn=0:1:N2-1;xn=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+sin(2*pi*f3*n/fs);%滿足條件的正弦序列Xk=fft(xn);AXk=abs(Xk(1:N2/2));%對x(n)進行128點DFT運算figure(4)k=(0:N2/2-1)*fs/N2;subplot(2,1,2)plot(k,AXk);%繪出X(k)的128點DFT頻譜幅度結(jié)果與仿真波形圖如下：圖①64點采樣正弦序列x(n)圖②X(k)64點DFT圖③X(k)補零到128點DFT圖④X(k)128點DFT因為采樣最小點數(shù)為82點，所以當(dāng)N1=64點采樣時進過D