数字信号处理实验报告

北京信息科技大学实验报告课程名称：数字信号处理实验项目：IIR数字滤波器设计实验仪器：计算机，MATLAB软件学院：信息与通信工程学院班级：姓名：学号：实验日期：指导老师：实验成绩：

实验二IIR数字滤波器设计一、实验目的1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。2、掌握数字滤波器的计算机软件实现方法。3、通过观察对实际心电图信号的滤波作用，学习数字滤波器在实际中的应用。二、实验仪器及材料计算机，MATLAB软件三、实验内容及要求1．设计巴特沃斯低通数字滤波器对人体心电信号进行滤波(1)人体心电图信号在测量过程中会受到工业高频干扰，所以必须经过低通滤波处理，才能作为判断心脏功能的有用信息。以下为一个实际心电图信号采样序列x(n)，其中存在高频干扰，采样周期T=1秒。在实验中，以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)用冲激响应不变法设计一个巴特沃斯低通IIR数字滤波器H(z)。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π,π]频率区间上，最小衰减大于15dB。写出数字滤波器H(z)的表达式，画出滤波器的幅频响应曲线。(3)用所设计的滤波器对实际心电图信号采样序列进行滤波处理，编写程序,计算H(z)对心电图信号采样序列x(n)滤波后的序列y(n)，并分别画出滤波前后的心电图信号波形图和频谱图。2．用help查看内部函数cheb1ord.m及cheby1.m，了解调用格式，并用此函数设计一个数字切贝雪夫带通滤波器。设计指标参数为：抽样频率fs=2kHz；在200Hz≤f≤400Hz时，最大衰减小于2dB；在f≤100Hz，f≥600Hz，最小衰减大于40dB。编程设计，求数字滤波器H(z)的表达式，画出滤波器的幅频响应曲线。四、实验程序(改为画五个图的程序，加粗的是查的函数的注释)%................实验二IIR数字滤波器设计....................%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)%画第一个图%FIGURE,byitself,createsanewfigurewindow,andreturnsitshandle.m=0:1:length(x)-1;%0到x(n)序列长度N-1的N个点，每隔1画一个点%length(x):求序列x(n)的长度subplot(1,2,1);%SUBPLOT(m,n,p,PROP1,VALUE1,…)setsthespecifiedpropertyvaluepairsonthe%subplotaxis.Toaddthesubplottoaspecificfigurepassthefigurehandleas%thevalueforthe'Parent'property.画图，是1行2列中的第一个图stem(m,x,'r');%以m为横坐标，x为纵坐标画离散曲线，曲线颜色为红色xlabel('时间序列n');ylabel('x(n)');%xlabel给x轴加标记“时间序列n”,ylabel给y轴加标记“x(n)”LL=0:27;X=fft(x,56);%将序列x(n)做56点快速傅立叶变换，结果赋给序列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横坐标范围是0-28，纵坐标范围是0-450wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;%先设计模拟滤波器，再冲激响应不变法变换为数字滤波器OmigaP=wp/T;OmigaS=ws/T;[N,Wn]=buttord(OmigaP,OmigaS,Rp,As,'s');%计算模拟butterworth滤波器的最小阶次n和截止频率为Wn.[numS,denS]=butter(N,Wn,'s');%设计模拟截止频率为Wn的n阶Butterworth低通滤波器，返回值为数字滤波器的系数。[numZ,denZ]=impinvar(numS,denS,1/T),%冲激响应不变法，返回值为数字滤波器的系数。%显示数字滤波器阶次与截止频率参数disp('模拟巴特沃斯滤波器的阶数N=');disp(N);disp('模拟巴特沃斯滤波器的3dB截止频率Wn=');disp(Wn);%DISP(X)displaysthearray,withoutprintingthearrayname，显示y=filter(numZ,denZ,x);%Y=FILTER(B,A,X)filtersthedatainvectorXwiththefilterdescribedbyvectors%AandBtocreatethefiltereddataY.%是设计出系数为B和A的数字滤波器对序列x（n）进行处理figure(2)%画第二个图subplot(1,2,1)stem(m,y,'.');xlabel('时间序列n');ylabel('y(n)')Y=abs(fft(y,56));%对序列fft(y,56)每一项取绝对值或模值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);%求频率响应函数，返回值为滤波器的频率响应db=20*log10((abs(Hz)+eps)/max(abs(Hz)));%max求最大值函数figure(3)%画第三个图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.开根号函数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('数字切贝雪夫带通滤波器的阶数N=');disp(2*n);%disp('数字切贝雪夫带通滤波器的通带截止频率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;%画第四个图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;%画第五个图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)五、实验结果图1图2.左为序列x(n)的时域图，右为序列x(n)的56点快速傅里叶变换（滤波前）图3.左为序列x(n)的时域图，右为序列x(n)的56点快速傅里叶变换（滤波后）图4.左为巴特沃兹低通滤波器在0-π之间的频率响应，右为以dB表示的频率响应图5.左为切比雪肤带通滤波器的频率响应，右为以dB表示的频率响应图6.左为巴特沃兹低通滤波器的零极点图