數(shù)字信號處理實(shí)驗(yàn)報(bào)告

北京信息科技大學(xué)實(shí)驗(yàn)報(bào)告課程名稱：數(shù)字信號處理實(shí)驗(yàn)項(xiàng)目：IIR數(shù)字濾波器設(shè)計(jì)實(shí)驗(yàn)儀器：計(jì)算機(jī)，MATLAB軟件學(xué)院：信息與通信工程學(xué)院班級：姓名：學(xué)號：實(shí)驗(yàn)日期：指導(dǎo)老師：實(shí)驗(yàn)成績：

實(shí)驗(yàn)二IIR數(shù)字濾波器設(shè)計(jì)一、實(shí)驗(yàn)?zāi)康?、熟悉用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的原理與方法。2、掌握數(shù)字濾波器的計(jì)算機(jī)軟件實(shí)現(xiàn)方法。3、通過觀察對實(shí)際心電圖信號的濾波作用，學(xué)習(xí)數(shù)字濾波器在實(shí)際中的應(yīng)用。二、實(shí)驗(yàn)儀器及材料計(jì)算機(jī)，MATLAB軟件三、實(shí)驗(yàn)內(nèi)容及要求1．設(shè)計(jì)巴特沃斯低通數(shù)字濾波器對人體心電信號進(jìn)行濾波(1)人體心電圖信號在測量過程中會(huì)受到工業(yè)高頻干擾，所以必須經(jīng)過低通濾波處理，才能作為判斷心臟功能的有用信息。以下為一個(gè)實(shí)際心電圖信號采樣序列x(n)，其中存在高頻干擾，采樣周期T=1秒。在實(shí)驗(yàn)中，以x(n)作為輸入序列，濾除其中干擾成分。x(n)=[-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]對序列x(n)用FFT做頻譜分析，生成x(n)的頻譜圖。(2)用沖激響應(yīng)不變法設(shè)計(jì)一個(gè)巴特沃斯低通IIR數(shù)字濾波器H(z)。設(shè)計(jì)指標(biāo)參數(shù)為：在通帶內(nèi)頻率低于0.2π時(shí)，最大衰減小于1dB；在阻帶內(nèi)[0.3π,π]頻率區(qū)間上，最小衰減大于15dB。寫出數(shù)字濾波器H(z)的表達(dá)式，畫出濾波器的幅頻響應(yīng)曲線。(3)用所設(shè)計(jì)的濾波器對實(shí)際心電圖信號采樣序列進(jìn)行濾波處理，編寫程序,計(jì)算H(z)對心電圖信號采樣序列x(n)濾波后的序列y(n)，并分別畫出濾波前后的心電圖信號波形圖和頻譜圖。2．用help查看內(nèi)部函數(shù)cheb1ord.m及cheby1.m，了解調(diào)用格式，并用此函數(shù)設(shè)計(jì)一個(gè)數(shù)字切貝雪夫帶通濾波器。設(shè)計(jì)指標(biāo)參數(shù)為：抽樣頻率fs=2kHz；在200Hz≤f≤400Hz時(shí)，最大衰減小于2dB；在f≤100Hz，f≥600Hz，最小衰減大于40dB。編程設(shè)計(jì)，求數(shù)字濾波器H(z)的表達(dá)式，畫出濾波器的幅頻響應(yīng)曲線。四、實(shí)驗(yàn)程序(改為畫五個(gè)圖的程序，加粗的是查的函數(shù)的注釋)%................實(shí)驗(yàn)二IIR數(shù)字濾波器設(shè)計(jì)....................%clc,clear,closeallT=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];figure(1)%畫第一個(gè)圖%FIGURE,byitself,createsanewfigurewindow,andreturnsitshandle.m=0:1:length(x)-1;%0到x(n)序列長度N-1的N個(gè)點(diǎn)，每隔1畫一個(gè)點(diǎn)%length(x):求序列x(n)的長度subplot(1,2,1);%SUBPLOT(m,n,p,PROP1,VALUE1,…)setsthespecifiedpropertyvaluepairsonthe%subplotaxis.Toaddthesubplottoaspecificfigurepassthefigurehandleas%thevalueforthe'Parent'property.畫圖，是1行2列中的第一個(gè)圖stem(m,x,'r');%以m為橫坐標(biāo)，x為縱坐標(biāo)畫離散曲線，曲線顏色為紅色xlabel('時(shí)間序列n');ylabel('x(n)');%xlabel給x軸加標(biāo)記“時(shí)間序列n”,ylabel給y軸加標(biāo)記“x(n)”LL=0:27;X=fft(x,56);%將序列x(n)做56點(diǎn)快速傅立葉變換，結(jié)果賦給序列y(n)XX=[X(1:28)];subplot(1,2,2);stem(LL,abs(XX),'.');xlabel('k=0,1,2,...,27');ylabel('|X(k)|');axis([0,28,0,450]);%axis橫坐標(biāo)范圍是0-28，縱坐標(biāo)范圍是0-450wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;%先設(shè)計(jì)模擬濾波器，再?zèng)_激響應(yīng)不變法變換為數(shù)字濾波器OmigaP=wp/T;OmigaS=ws/T;[N,Wn]=buttord(OmigaP,OmigaS,Rp,As,'s');%計(jì)算模擬butterworth濾波器的最小階次n和截止頻率為Wn.[numS,denS]=butter(N,Wn,'s');%設(shè)計(jì)模擬截止頻率為Wn的n階Butterworth低通濾波器，返回值為數(shù)字濾波器的系數(shù)。[numZ,denZ]=impinvar(numS,denS,1/T),%沖激響應(yīng)不變法，返回值為數(shù)字濾波器的系數(shù)。%顯示數(shù)字濾波器階次與截止頻率參數(shù)disp('模擬巴特沃斯濾波器的階數(shù)N=');disp(N);disp('模擬巴特沃斯濾波器的3dB截止頻率Wn=');disp(Wn);%DISP(X)displaysthearray,withoutprintingthearrayname，顯示y=filter(numZ,denZ,x);%Y=FILTER(B,A,X)filtersthedatainvectorXwiththefilterdescribedbyvectors%AandBtocreatethefiltereddataY.%是設(shè)計(jì)出系數(shù)為B和A的數(shù)字濾波器對序列x（n）進(jìn)行處理figure(2)%畫第二個(gè)圖subplot(1,2,1)stem(m,y,'.');xlabel('時(shí)間序列n');ylabel('y(n)')Y=abs(fft(y,56));%對序列fft(y,56)每一項(xiàng)取絕對值或模值subplot(1,2,2)stem(LL,abs(Y(1:28)),'.');xlabel('k=0,1,2,...,27');ylabel('|Y(k)|');axis([0,28,0,450]);w=0:0.01:3.14;Hz=freqz(numZ,denZ,w);%求頻率響應(yīng)函數(shù)，返回值為濾波器的頻率響應(yīng)db=20*log10((abs(Hz)+eps)/max(abs(Hz)));%max求最大值函數(shù)figure(3)%畫第三個(gè)圖subplot(1,2,1),plot(w/pi,abs(Hz),'r');xlabel('frequencyinpiunits');ylabel('|H(e(jw))|');axis([0,1,0,1.1]);subplot(1,2,2),plot(w/pi,db,'r');xlabel('frequencyinpiunits');ylabel('Decibels');axis([0,1,-40,2]);set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);set(gca,'YTickMode','manual','YTick',[-As,-Rp,0]);grid;%SET(H,'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2,...)%setsmultiplepropertyvalueswithasinglestatement.Notethatitispermissible%touseproperty/valuestringpairs,structures,andproperty/valuecellarray%pairsinthesamecalltoSET.%GRIDGridlines.fs=2000;Rp=2;wp1=2*pi*200/fs;wp2=2*pi*400/fs;Wp=[wp1/pi,wp2/pi];As=40;ws1=2*pi*100/fs;ws2=2*pi*600/fs;Ws=[ws1/pi,ws2/pi];ep=sqrt(10^(Rp/10)-1);%SQRT(X)isthesquarerootoftheelementsofX.開根號函數(shù)Ripple=sqrt(1/(1+ep*ep));Attn=1/(10^(As/20));disp('**********************************************************')[n,Wn]=cheb1ord(Wp,Ws,Rp,As);%[N,Wn]=CHEB1ORD(Wp,Ws,Rp,Rs)returnstheorderNofthelowest%orderdigitalChebyshevTypeIfilterthatlosesnomorethanRpdB%inthepassbandandhasatleastRsdBofattenuationinthestopband.%WpandWsarethepassbandandstopbandedgefrequencies,normalized%from0to1(where1correspondstopiradians/sample).%disp('數(shù)字切貝雪夫帶通濾波器的階數(shù)N=');disp(2*n);%disp('數(shù)字切貝雪夫帶通濾波器的通帶截止頻率Wn=');disp(Wn);[b,a]=cheby1(n,Rp,Wn,'bandpass')%[B,A]=CHEBY1(N,R,Wp)designsanNthorderlowpassdigitalChebyshevfilterwith%Rdecibelsofpeak-to-peakrippleinthepassband.CHEBY1returnsthefilter%coefficientsinlengthN+1vectorsB(numerator)andA(denominator).%Thepassband-edgefrequencyWpmustbe0.0<Wp<1.0,with1.0correspondingto%halfthesamplerate.[Hw]=freqz(b,a);%[H,W]=FREQZ(B,A,N)returnstheN-pointcomplexfrequencyresponsevectorHandtheN-pointfrequencyvectorWinradians/sampleofthefilter:jw-jw-jmwjwB(e)b(1)+b(2)e+....+b(m+1)eH(e)=----=------------------------------------jw-jw-jnwA(e)a(1)+a(2)e+....+a(n+1)egivennumeratoranddenominatorcoefficientsinvectorsBandA.ThefrequencyresponseisevaluatedatNpointsequallyspacedaroundtheupperhalfoftheunitcircle.IfNisn'tspecified,itdefaultsto512.db=20*log10((abs(H)+eps)/max(abs(H)));%plotfigure(4);clf;%畫第四個(gè)圖subplot(1,2,1);plot(w/pi,abs(H));title('BandpassFilterMagnitudeResponse');xlabel('frequencyinpiunits');ylabel('|H(e(jw))|');axis([0,1,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,0.10.20.30.40.61]);set(gca,'YTickMode','manual','YTick',[0,Attn,Ripple,1]);grid;subplot(1,2,2);plot(w/pi,db);title('BandpassFilterMagnitudeindB');xlabel('frequencyinpiunits');ylabel('Decibels');axis([0,1,-60,2]);set(gca,'XTickMode','manual','XTick',[0,0.10.20.30.40.61]);set(gca,'YTickMode','manual','YTick',[-As,-Rp,0]);grid;figure(5);clf;%畫第五個(gè)圖hf=findobj(gcf,'Color','b');set(hf,'Color','r','linewidth',1.5);%H=FINDOBJ('P1Name',P1Value,...)returnsthehandlesoftheobjectsattheroot%levelandbelowwhosepropertyvaluesmatchthosepassedasparam-valuepairs%totheFINDOBcommand.subplot(1,2,1);zplane(b,a);title('Zeros&Poles');subplot(1,2,2);[ht]=impz(b,a);h=(h([1:1:150]));stem(h,'.');title('ImpulseResponse');axis([0150-0.20.2])hf=findobj(gcf,'Color','b');set(hf,'linewidth',1.5)五、實(shí)驗(yàn)結(jié)果圖1圖2.左為序列x(n)的時(shí)域圖，右為序列x(n)的56點(diǎn)快速傅里葉變換（濾波前）圖3.左為序列x(n)的時(shí)域圖，右為序列x(n)的56點(diǎn)快速傅里葉變換（濾波后）圖4.左為巴特沃茲低通濾波器在0-π之間的頻率響應(yīng)，右為以dB表示的頻率響應(yīng)圖5.左為切比雪膚帶通濾波器的頻率響應(yīng)，右為以dB表示的頻率響應(yīng)圖6.左為巴特沃茲低通濾波器的零極點(diǎn)圖