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

测试信号分析与处理实 验 指 导 书翟任何 王明武 编 写适用专业：测控专业陕西理工学院机械工程学院二零壹零年十月目 录实验一：数字滤波器的设计第2页实验二：数字振荡器的实现第14页实验一：数字滤波器的设计实验学时：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平面上，得到：（2-3）（二）用窗函数设计FIR数字低通滤波器如果所希望的滤波器的理想频率响应函数为 ，则其对应的单位脉冲响应为 = （2-4）窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，且是非因果的，所以用窗函数将截断，并进行加权处理，得到： (2-5)作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为： (2-6)式中，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是否满足要求，要进行验算。一般在尾部加零使长度满足于2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足 (2-4)根据上式中的正负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据设计的滤波特性正确选择其中一类。四、实验条件：微型计算机 matlab软件五、实验步骤：1、复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容。根据设计指标，调用MATLAB信号处理工具箱函数buttord和butter，可得到H（z）。编写滤波器仿真程序，计算H(z)对心电图信号采样序列x(n)的相应序列y(n)。在通过计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验内容。2、复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理掌握设计步骤。编写能产生矩型窗、哈明窗、汉宁窗、莱克曼窗的窗函数子程序。 编写主程序。主程序框图如图21所示，仅供参考。其中幅度特性要求用分贝dB表示。开始读入窗口长度N计算hd（n）调窗函数子程序求（n）计算h（n）=hd（n）*（n）调FFT子程序对h（n）进行DFT调绘图子程序绘制幅频特性曲线相位特性曲线结束计算幅度特性和相位特性图21 主程序框图设 ： (2-4) (2-5) (2-6)画图时，20lg打印幅度特性。第k点对应的频率。为使曲线包络更接近的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在末尾补零的方法，使长度变为64，再进行4点DFT，则可得到更精确的幅度衰减特性曲线。 六、思考题：1、用双线性变换法设计数字滤波器过程中,变换公式：s=中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, 0 0, 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;stem(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)后的对数频率响应20logHa3(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(20logHa3(ejw);title(通过滤波器H1(z)、H2(z)、H3(z)后的对数频率响应20logHa3(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;h=sin(w*m)./(pi*m);if k=1 B=bartlett(N);else if k=2 B=hanning(N); else if k=3 B=hamming(N); else if k=4 B=blackman(N); end end endendhd=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(20logH(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;close all;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; N=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); if k=1 B=boxcar(N); string=Boxcar,N=,num2str(N); else if k=2 B=hamming(N); string=Hamming,N=,num2str(N); else if k=3 B=hanning(N); string=Hanning,N=,num2str(N); else if k=4 B=blackman(N); string=Blackman,N=,num2str(N); end end end end h=hd.*(B); %得到FIR数字滤波器 H,m=freqz(h,1,1024,whole); %求其频率响应 mag=abs(H); %得到幅值 db=20*log10(mag+eps)/max(mag); pha=angle(H); %得到相位 i=i+1; figure(i) subplot(2,2,1); n=0:N-1; stem(n,h,.); axis(0,N-1,-0.1,0.3); hold on; n=0:N-1; x=zeros(N); plot(n,x,-); xlabel(n); ylabel(h(n); title(实际低通滤波器的h(n); text(0.3*N),0.27,string); hold off; subplot(2,2,2); plot(m/pi,db); axis(0,1,-100,0); xlabel(w/pi); ylabel(dB); title(衰减特性(dB); grid; subplot(2,2,3); plot(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,-); title(相频特性); xlabel(频率(rad); ylabel(相位(rad); axis(0,3.15,-4,4); subplot(2,2,4); plot(m,mag); title(频率特性); xlabel(频率W(rad); ylabel(幅值); axis(0,3.15,0,1.5); text(0.9,1.2,string); b=menu(Do You want To Continue ?,Yes,No); if b=2 b=0; endendtemp=menu(Close All Figure ?,Yes,No);if temp=1 close allend实验二：数字振荡器的实现实验学时：2实验类型：常规实验要求：必修一、实验目的：1、 学习数字振荡器的DSP实现原理2、 学习C54X定时器的使用3、 学习C54X中断程序的编写二、实验内容：利用定时器产生一个2kHz的正弦信号。使用汇编语言和C语言分别完成源程序的编写，并利用CCS的图形显示功能查看波形。三、实验原理：1、 正弦余弦信号的产生方法：在数字信号处理中，会经常使用到正弦余弦信号。通常的方法是将某个频率的余弦值预先计算出来后制成一个表，DSP工作时仅作查表运算即可。另外一种获得正弦余弦信号的方法是：利用数字振荡器用叠代方法产生正弦信号。定时器设置成每25uS产生一次中断（等效于采样速率为40K）。利用该中断，在中断服务程序中用叠代算法计算出一个SIN值。2、 程序实现：(1) 数字振荡器原理：设一个传递函数为正弦序列sinkT，其Z变换为：H(Z)=其中，A=2cosT,B=-1,C=sinT。设初始条件为0，求出上式的反变换位得：ykAyk-1Byk-2Cxk-1这是一个二阶差分方程，其单位冲击相应即为sinkT。利用单位冲击函数xk-1的性质，即仅当时k=1时，xk-1=1，代入上式得：k=0 y0Ay-1By-20 = 0k=1 y1Ay1By-1c = Ck=2 y2Ay1By00 = ACk=3 y3Ay2By1k=n ynAyn-1Byn-2在k2以后，yk能用yk-1和yk-2算出，这是一个递归的差分方程。(2) C54X的定时器操作C54X的片内定时器利用CLKOUT时钟计数，用户使用三个寄存器（TIM，PRD，TCR）来控制定时器，表中列出了定时器控制各个寄存器的各比特位的具体定义。VC5402的另一个定时器（定时器1）的控制寄存器分别为：0x30（TIM1），0x31（PRD1），0x32（TCR1）。寄存器地址名称用途0024hTIM定时器寄存器，每计数一次，自动减10025hPRD定时器周期寄存器，当TIM减为0后，CPU自动将PRD的值装入TIM0026hTCR定时器控制寄存器(3) C54X中断的使用在C54X中用户可以通过中断屏蔽寄存器来决定开放或关闭一个中断请求。寄存器的各个比特位的定义如表：15、1413121110987ResvdDMAC5DMAC4BXINT1 orDMAC3BXINT1 orDMAC2HPINTINT3TINTorDMAC16543210DMAC0BXINT0TINT0TINT0INT2INT1INT0其中，HPINT表示HPI接口中断，INT3-INT0为外部引脚产生的中断，TXINT和TRINT为TDM串口的发送和接收中断，BXINT0和BRINT0为BSP串口的发送和接收中断，TINT0定时器0中断。在中断屏蔽寄存器IMR中，1表示允许CPU响应对应的中断，0表示禁止。当然要CPU响应中断，ST1寄存器中的INTM还应该为0（允许所有中断）。当DSP响应中断时，PC指针指向中断向量表中对应中断的地址，进入中断服务子程序。中断向量表是C54X存放中断服务程序的一段内存区域，大小为80H。在中断向量表中，每一个中断占用4个字的空间，一般情况是将一条跳转或延时跳转指令存放于此。当然，如果中断服务程序很短（小于或等于4个字），可以直接放入该向量表。中断向量表的位置可以通过修改基地址来改变，其基地址由PMST寄存器中的IPTR决定。四、实验条件：DSP实验箱、装有CCS环境的计算机五、实验步骤：a) 使用汇编语言编制程序1）启动CCS，新建工程文件sinewave.pjt。2）编写汇编源程序exer3.asm。3）选择Project 菜单中的Add File to Project 选项， 将汇编源程序exer3.asm、vec_table.asm 和连接定位sinewave.cmd 文件依次添加到工程文件中。其中，exer3.asm 包括初始化代码和中断服务程序，而vec_talbe.asm 包含中断向量表。注意：你可以在添加文件对话框中选择显示不同的文件类型来加快文件选择速度。4）选择Project 菜单中的Options 选项，或使用鼠标右键单击工程文件名（如sinewave.mak）并选择Options 项来修改或添加编译、连接中使用的参数。例如，选择Assembler 窗口，选择“Enable Symbolic Debug Infomation”以便使用汇编源代码级调试（你可以在汇编源程序设置断点等等）。选择Linker 窗口，在“Output Filename”栏中写入输出OUT 文件的名字，如sine.out，你还可以设置生成的MAP 文件名。5）完成编译、连接，正确生成OUT 文件。然后使用File 菜单的“Load Program”选项，将生成的OUT 文件（如sine.out）装入DSP 的片内存储器。这时CCS 将显示程序的启始地址_c_int00。6）选ViewGraphTime/Frequency打开图形显示设置窗口。在弹出的对话框中按下图设置，主要修改“Start Address”为y0（y0 为生成的正弦波输出变量）；“ Acquisition Buffer Size ” 为1 ， “ DSP Data Type ” 为“ 16-bit signedinteger”。7）在汇编源程序的中断服务程序（_tint）中的“nop”语句处设置断点。该行被加亮为洋红色。选择DebugAnimate，运行程序，观察输出波形。8）用右键单击图形显示窗口，并选择“Proporties”项以便修改显示属性。将“ Display Type ” 项改为“ FFT Magnitude ” 以便显示信号频谱。修改“ SamplingRate(Hz)”项为40000，然后退出。注意观察生成的正弦波频率。9）清除所有断点，关闭除波形显示窗口外的所有窗口，并关闭工程文件。b) 使用C语言编制程序1）新建工程文件sinewave_c.pjt。2）编写汇编源程序timer.c。3）添加timer.c、vec_table.asm 源程序，再添加timer.cmd，再添加C 使用的标准库rts.lib。该文件应该在CCS 安装目录中。例如，若CCS 安装在d:ti 下，则rts.lib 应该在d:tic5400cgtoolslib下。4）完成编译、连接，正确生成OUT 文件。3）将生成的OUT 文件（如timer.out）装入DSP 的片内存储器。这时CCS 将显示程序的启始地址_c_int00。4）打开C 源程序（timer.c）窗口，在中断服务程序（函数tint()）的“con_buf=0;”语句处增加一个断点。同样打开图形显示窗口，并将“Start Address”改为buf；“Acquisition Buffer Size”改为128，“Display Data Size”改为128，“DSP Data Type ” 为“ 32-bit floating point ”。5）选择DebugAnimate，运行程序，观察输出波形。同样用右键单击图形显示窗口，显示信号频谱。注意修改“Sampling Rate(Hz)”项为40000，然后退出。注意观察生成的正弦波频率。6）清除所有断点，关闭除波形显示窗口外的所有窗口，并关闭工程文件。六、思考题：4、 为什么在频谱图中会观察到高频成分？5、 改变数字振荡器参数，使采样频率为20KHZ，输出正弦信号频率为4KHZ七、实验报告：（一）标题栏（二）实验目的（三）实验内容（四）源代码（五）运行结果（六）问题回答（七）实验总结八、源代码：exer3.asm ：.title for test INT service program .(25 us) .mmregs .global _c_int00,_tint,vector OFF_INTIMER .set 04Ch ;INIT_A .set 079bch ; A/2=0.9510498INIT_B .set 0c000h ; B/2=-0.5INIT_C .set 013c7h ; C/2=0.1545105 .bss y0,1 .bss y1,1 .bss y2,1 .bss temp,1 .bss AA,1 .bss BB,1 .bss CC,1 .text _c_int00: ld #0,dp st #1fffh,sp ssbx intm ; disable all interrupt ! ld #vector, a ; get vector table address ! and #0FF80h, a andm #007Fh, pmst or pmst, a stlm a, pmst ; 设置IPTR stm #10h,TCR ; init TIMER stm #2499,PRD ; f=100M/(2499+1)=40kHz stm #20h,TCR ; reset TIMER ldm IMR,a ; read back IMR or #08h,a ; enable TIMER interrupt stlm a,IMR ; set IMR ld #temp,dp ; set DP ssbx FRCT ; prepare for fraction mpy st #INIT_A,AA ; init AA,BB,CC st #INIT_B,BB ; st #INIT_C,CC ; pshd CC popd y2 ; init y2,y2=CC ld AA,T ; T=AA mpy y2,a ; y2*AA - a sth a,y1 ; y2*AA - y1 stm #0h,TCR ; enable TIMER nop rsbx intm ; enable all int !again: nop b again nop nop nop nop nop nop;-; interrupt for INT_TIMER !;-_tint: ld #BB,DP ld BB,T ; T=BB mpy y2,a ; a=y2*BB ltd y1 ; T=y1,y2=y1