测试信号分析与处理实验指导书

《测试信号分析与处理》实验指导书翟任何王明武编写适用专业：测控专业陕西理工学院机械工程学院二零壹零年十月1目录实验一：数字滤波器的设计·······················第2页实验二：数字振荡器的实现······················第14页2实验一：数字滤波器的设计实验学时：4实验类型：综合实验要求：必修一、实验目的：1、熟悉Matlab界面并进行操作。2、掌握数字滤波器的计算机仿真方法。2、掌握用双线性变换法设计IIR数字滤波器的原理与方法，通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。3、掌握用窗函数法设计FIR数字滤波器的原理与方法，了解各种窗函数对滤波特性的影响。二、实验内容：（一）IIR数字滤波器设计（1）用双线性变换法设计一个巴特沃斯低通IIR数字滤波器，设计指标参数为：在通带内频率低于0.2pi时,最大衰减小于1dB;在阻带内[0.3pi,pi]频率区间上，,最小衰减大于15dB；（2）以0.02pi为采样间隔，打印出数字滤波器在频率区间[0,0.5pi]上的幅频响应特性曲线；（3）用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图波形图，观察总结滤波作用与效果。（二）FIR数字滤波器设计（1）用四种窗函数设计线性相位低通FIR数字滤波器，截止频率π/4rad，N=256。（2）绘制相应的幅频特性曲线，观察3dB带宽和20dB带宽以及阻带最小衰减。（3）比较四种窗函数对滤波器特性的影响。三、实验原理：（一）用双线形变换法设计IIR数字低通滤波器脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性，产生的原因是模拟低通滤波器不是带限于折叠频率π/T，在数字化后产生了频谱混叠，再通过映射关系，使数字滤波器在ω=π附近形成频谱混叠现象。为了克服这一缺点，可以采用非线性频率压缩方法，将整个模拟频率轴压缩到±π/T之间，再用转换到Z平面上。设，然后经过非线性频率轴压缩用表示。用正切变换实现频率压缩（2-1）式中T仍是采样间隔，当Ω1从π/T经过0变化到π/T时，Ω则由∞经过0变化到+∞，实现了s平面上整个虚轴完全压缩到s1平面上虚轴的±π/T之间的转换。这样便有（2-2）再通过转换到z平面上，得到：121tan()2TT1112121()21sTsTesthTTTe1sTzesTzesTze(),aHssj111ˆ(),aHssj3（2-3）（二）用窗函数设计FIR数字低通滤波器如果所希望的滤波器的理想频率响应函数为)(jdeH，则其对应的单位脉冲响应为)(nhd=21deeHjjd)(（2-4）窗函数设计法的基本原理是用有限长单位脉冲响应序列)(nh逼近)(nhd。由于)(nhd往往是无限长序列，且是非因果的，所以用窗函数)(n将)(nhd截断，并进行加权处理，得到：)(nh＝)(nhd)(n(2-5))(nh作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数)(jdeH为：)(jdeH＝10)(Nnjenh(2-6)式中，N为所选窗函数)(n的长度。窗函数法设计的滤波器性能取决于窗函数)(n的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。表2.1各种窗函数的基本参数窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134π/N-12三角形窗-258π/N-25汉宁窗-318π/N-44哈明窗-418π/N-53不莱克曼窗-5712π/N-74凯塞窗(α=7.865)-5710π/N-80)(jeH是否满足要求，要进行验算。一般在)(nh尾部加零使长度满足于2的整数次幂，以便用FFT计算)(jeH。如果要观察细节，补零点数增多即可。如果)(jeH不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。如果要求线性相位特性，则)(nh还必须满足)1()(nNhnh(2-4)根据上式中的正负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据设计的22sTzsT4滤波特性正确选择其中一类。四、实验条件：微型计算机matlab软件五、实验步骤：1、复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容。根据设计指标，调用MATLAB信号处理工具箱函数buttord和butter，可得到H（z）。编写滤波器仿真程序，计算H(z)对心电图信号采样序列x(n)的相应序列y(n)。在通过计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验内容。2、复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理掌握设计步骤。①编写能产生矩型窗、哈明窗、汉宁窗、莱克曼窗的窗函数子程序。②编写主程序。主程序框图如图2－1所示，仅供参考。其中幅度特性要求用分贝dB表示。图2－1主程序框图开始读入窗口长度N计算hd（n）调窗函数子程序求ω（n）计算h（n）=hd（n）*ω（n）调FFT子程序对h（n）进行DFT调绘图子程序绘制幅频特性曲线相位特性曲线结束计算幅度特性和相位特性5设：)()(nhDFTkH(2-4))()()(KjHkHkHIR(2-5))()()(22kHkHkHIR(2-6)画图时，20lg)(kH打印幅度特性。第k点对应的频率kNk2。为使曲线包络更接近)(jeH的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在)(nh末尾补零的方法，使长度变为64，再进行4点DFT，则可得到更精确的幅度衰减特性曲线。六、思考题：1、用双线性变换法设计数字滤波器过程中,变换公式：s=11z1z1T2中T的取值,对设计结果有无影响?为什么?2、如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出设计步骤。七、实验报告：（一）标题栏（二）实验目的（三）实验内容（四）源代码（五）运行结果（六）问题回答（七）实验总结八、源代码：参考程序1：人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样式本x(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,00,0,-2,-4,0,0,0,-2,-2,0，0,-2,,-2,-2,-2,0}参考程序：%x(n)的心电脉冲函数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];subplot(2,2,1);n=0:55;6stem(n,x,'.');xlabel('n');ylabel('x(n)');title('x(n)的心电脉冲函数');%通过滤波器H1(z)后的y1(n)函数A=0.09036;b1=[A,2*A,A];a1=[1,-1.2686,0.7051];h1=filter(b1,a1,x);[H1,w]=freqz(b1,a1,100);%通过滤波器H1(z)、H2(z)后的y2(n)函数b2=[A,2*A,A];a2=[1,-1.0106,0.3583];h2=filter(b2,a2,h1);[H2,w]=freqz(b2,a2,100);%通过滤波器H1(z)、H2(z)、H3(z)后的y3(n)函数b3=[A,2*A,A];a3=[1,-0.9044,0.2155];h3=filter(b3,a3,h2);[H3,w]=freqz(b3,a3,100);subplot(2,2,2);stem(n,h3,'.');xlabel('n');ylabel('y(n)');title('通过滤波器H1(z)、H2(z)、H3(z)后的y3(n)函数');%通过滤波器H1(z)、H2(z)、H3(z)后的对数频率响应20log[Ha3(ejw)]函数subplot(2,2,3);H4=H1.*(H2);H=H4.*(H3);mag=abs(H);db=20*log10((mag+eps)/max(mag));plot(w/pi,db);xlabel('w/п');ylabel('20log[Ha3(ejw)]');title('通过滤波器H1(z)、H2(z)、H3(z)后的对数频率响应20log[Ha3(ejw)]函数');grid;参考程序2：N=input('窗宽度N=');k=input('窗型：1.矩形窗,2.hanning(升余弦窗),3.hamming（改进的升余弦窗）,4.Blackman请选择：');subplot(2,2,1);w=pi/5;a=(N-1)/2;n=0:(N-1);m=n-a+eps;7h=sin(w*m)./(pi*m);ifk==1B=bartlett(N);elseifk==2B=hanning(N);elseifk==3B=hamming(N);elseifk==4B=blackman(N);endendendendhd=h.*(B');stem(n,hd,'.');xlabel('n');ylabel('h(n)');title('在矩形窗下的N=33时h(n)函数');subplot(2,2,2);[H,m]=freqz(hd,[1],1024);mag=abs(H);db=20*log10((mag+eps)/max(mag));plot(m/pi,db);xlabel('w/п');ylabel('20log[H(ejw)]');title('h(n)的幅频特性');grid;pha=angle(H);subplot(2,2,3);plot(m,pha);xlabel('n');ylabel('φ');title('h(n)的相频特性');subplot(2,2,4);plot(m,mag);xlabel('w');ylabel('H(ejw)');title('h(n)的幅频特性');参考程序3：b=1;closeall;i=0;while(b);temp=menu('选择窗函数长度N','N=10','N=15','N=20','N=25','N=30','N=33','N=35','N=40','N=45','N=50','N=55','N=60','N=64');menu1=[10,15,20,25,30,33,35,40,45,50,55,60,64];8N=menu1(temp);temp=menu('选择逼近理想低通滤波器截止频率Wc','Wc=pi/4','Wc=pi/2','Wc=3*pi/4','Wc=pi','Wc=0.5','Wc=1.0','Wc=1.5','Wc=2.0','Wc=2.5','Wc=3.0');menu2=[pi/4,pi/2,3*pi/4,pi,0.5,1,1.5,2,2.5,3];w=menu2(temp);n=[0:(N-1)];hd=ideal(w,N);%得到理想低通滤波器k=menu('请选择窗口类型:','boxcar','hamming','hanning','blackman');ifk==1B=boxcar(N);string=['Boxcar','N=',num2str(N)];elseifk==2B=hamming(N)