《数字信号处理》实验指导书

《数字信号处理》上机实验指导书 -《数字信号处理》上机实验指导书 --《数字信号处理》上机实验指导书一、引言“数字信号处理”是一门理论和实验密切结合的课程，为了深入地掌握课程内容，应当在学习理论的同时，做习题和上机实验。上机实验不仅可以帮助学生深入地理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。所以，根据本课程的重点要求编写了四个实验。第一章、二章是全书的基础内容，抽样定理、时域离散系统的时域和频域分析以及系统对输入信号的响应是重要的基本内容。由于第一、二章大部分内容已经在前期《信号与系统》课程中学习完，所以可通过实验一帮助学生温习以上重要内容，加深学生对“数字信号处理是通过对输入信号的一种运算达到处理目的” 这一重要概念的理解。这样便可以使学生从《信号与系统》课程顺利的过渡到本课程的学习上来。第三章、四章DFT、FFT是数字信号处理的重要数学工具，它有广泛的使用内容。限于实验课时，仅采用实验二“用FFT对信号进行谱分析”这一实验。通过该实验加深理解DFT的基本概念、基本性质。FFT是它的快速算法，必须学会使用。所以，学习完第三、四章后，可安排进行实验二。数字滤波器的基本理论和设计方法是数字信号处理技术的重要内容。学习这一部分时，应重点掌握IIR和FIR两种不同的数字滤波器的基本设计方法。IIR滤波器的单位冲激响应是无限长的，设计方法是先设计模拟滤波器，然后再通过 S~Z平面转换，求出相应的数字滤波器的系统函数。这里的平面转换有两种方法，即冲激响应不变法和双线性变换法，后者没有频率混叠的缺点，且转换简单，是一种普遍应用的方法。学习完第六章以后可以进行实验三。FIR滤波器的单位冲激响应是有限长的，设计滤波器的目的即是求出符合要求的单位冲激响应。窗函数法是一种基本的，也是一种重要的设计方法。学习完第七章后可以进行实验四。以上所提到的四个实验，可根据实验课时的多少恰当安排。例如：实验一可根据学生在学习《信号与系统》课程后，掌握的程度来确定是否做此实验。若时间紧，可以在实验三、四之中任做一个实验。建议实验二必做，因为DFT和FFT是学习数字信号处理必须掌握的内容。二、关于使用计算机语言由于数字信号处理实验的主要目的是验证数字信号处理的有关理论，进一步理解巩固所学理论知识。所以，实现实验用的算法语言可以有许多种，但为了提高实验效率，要求学生用编程效率比C语言高好几倍的MATLAB语言。下面介绍MATLAB的主要特点。（有关MATLAB的启动、程序运行和有关信号处理工具箱函数等内容将放到最后附录中介绍。）MATLAB是一种交互式的以矩阵为基本数据结构的系统。在生成矩阵对象时，不要求明确的维数说明。所谓交互式，是指MATLAB的草稿纸编程环境。即用户每输入一条命令并按回车键，MATLAB系统便解释执行之，并显示执行结果。根据该结果，用户立即知道刚输入的命令的正确性，或利用中间结果进行其他处理等。与C语言或FORTRON语言做科学数值计算的程序设计相比较，利用MATLAB可节省大量的编程时间。将其用于数字信号处理实验，则可大大提高实验效率，在有限的上机时间内，实验内容可增加几倍。例如， C语言FFT子程序有70多行，而用MATLAB只调用一个fft函数即可实现对序列进行FFT计算。另外，MATLAB的工具箱及图形显示（打印）功能，可满足各层次人员直观、方

便的进行分析、计算和设计工作，从而可大大节省时间。例如，序列的卷积、滤波，系统函数H(z)的幅频特性和相频特性等计算，均有现成的工具箱函数。而用其它算法语言完成这些计算的编程比较麻烦，且程序较长。由于上述特点，在美国一些大学里， MATLAB已成为辅助教学的有益工具。MATLAB已成功地用于数字信号处理课程中的问题分析、实验、滤波器设计及计算机模拟。附录中所介绍的信号处理工具箱函数及绘图函数基本可满足本教材所要求的上机实验需要。对序列进行谱分析的MATLAB程序及运行结果见附录。三、上机实验实验一：信号、系统及系统响应1.实验目的①熟悉连续信号经过理想抽样前后的频谱变化关系，加深对时域抽样定理的理解。②熟悉时域离散系统的时域特性。③利用卷积方法观察分析系统的时域特性。④掌握序列傅里叶变换的计算机实验方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。2.实验原理与方法抽样是连续信号数字处理的第一个关键环节。对抽样过程的研究不仅可以了解抽样前后信号时域和频域特性发生的变化以及信号信息不丢失的条件，而且可以加深对傅里叶变换、Z变换和序列傅里叶变换之间关系式的理解。我们知道，对一个连续信号Xa⑴进行理想抽样的过程可用(1.1)式表示。XaXa(t)=Xa(t)督t)(1.1)其中？a⑴为Xa(t)的理想抽样，标⑴为周期冲激脉冲，即T(t)(tnT)—— (1.2)nXa(t)的傅里叶变换兄(j)为)?a(j)=1Xa[j(ks)]——(1.3)

Tk(1.3)式表明)?a(j)为Xa(j)的周期延拓，其延拓周期为抽样角频率(Qs=2订)。抽样前后信号的频谱示意图见“参考教材数字信号处理教程(程佩青)图1-29”。只有满足抽样定理时，才不会发生频率混叠失真。数字信号处理教程(程佩青)在计算机上用高级语言编程直接按(1.3)式计算理想抽样？a(t)的频谱£a(j)很不方便。下面导出用序列的傅里叶变换来计算 £a(j)的公式将(1.2)式代入(1.1)式并进行傅里叶变换，&(j)= XUt)ejtdt= Xa(nT)(tnT)ejtdtnxa(nT)exa(nT)ejnT(tnT)dtnXa(nT)ejnTn(1.4)式中的xa(nT)就是采样后得到的序列x(n),即x(n)=Xa(nT)x(n)的序列傅里叶变换为X(j)= x(n)ejn——(1.5)n比较(1.5)和(1.4)可知父a(j)=X(ej)|a=QT (1.6)这说明两者之间只在频率度量上差一个常数因子 To实验过程中应注意这一差别。离散信号和系统在时域均可用序列来表示。序列图形给人以形象直观的印象，它可加深我们对信号和系统的时域特征的理解。本实验还将观察分析几种信号及系统的时域特性。为了在数字计算机上观察分析各种序列的频域特性，通常对 X(ejw)在［0,2耳上进行M点采样来观察分析。对长度为N的有限长序列x(n),有N1X(ejk) x(n)ejkn——(1.7)n0… 2.其中 k一k,k=0,1，…，M-1M通常M应取得大一些，以便观察谱的细节变化。取模|X(ejk)|可绘出幅频待性曲线。一个时域离散线性非移变系统的输入/输出关系为y(n)=x(n)*h(n)=x(m)h(nm) (1.8)m这里，y(n)为系统的输出序列，x(n)为输入序列。h(n)、x(n)可以是无限长，也可以是有限长。为了计算机绘图观察方便，主要讨论有限长情况。如果h(n)和x(n)的长度分别为N和M,则y(n)的长度为L=N+M-1。这样，(1.8)式所描述的卷积运算就是序列移位、相乘和累加的过程，所以编程十分简单。上述卷积运算也可以在频域实现(即卷积定理：时域卷积，频域相乘。)Y(ej)=X(ej，H(ej)——(1⑼(1.9)式右边的相乘是在各频点｛wk｝上的频谱值相乘。.实验内容(1)连续信号分析①连续时间非周期信号的选择参考：信号与系统(第二版)上册；郑君里一附录三 (P378)1单边指数脉冲；2双边指数脉冲；3钟形脉冲；4抽样脉冲；

②用MATLAB编制程序广生信号Xa⑴求连续时间非周期信号的傅氏变换Xa(jQ)1分别绘制Xa(t),|Xa(jf)|图形/1、C结束)

图1.1连续信号分析程序框图③时域观察，频域分析调整信号Xa(t)参数，观察时域波形变化。分析频域波形，选取近似带限频率值fmax,确定对Xa⑴—X(n)=Xa(nT)的抽样周期T(T=1/fs),不失真的条件应该满足抽样频率fs>2fmaX(即：时域抽样定理)。例如图1.2给出了单边减幅余弦信号Xa(t)=eXp(-10*t).*cos(2*pi*25*t)的时域和幅频特性曲线，由此图可以近似选取Xa⑴的抽样频率fs>2fmax=2x150=300Hz。图1.2Xa⑴的时域和幅频特性曲线(2)离散信号分析①离散时间非周期信号x(n)的生成对前一个实验中的连续信号Xa⑴抽样，用于产生实验中要用到的信号序列x(n)=Xa(nT)。根据上一个实验分析结果选取近似带限频率值 fmax,按照抽样频率fs=2fmax、fs>2fmax和fs<2fmax三种情况编制MATLAB源程序(T=1/fs)。

②用MATLAB编制程序图1.3离散信号分析程序框图③时域观察，频域分析调整抽样频率fs=2fmax、fs>2fmax和fs<2fmax,观察时域波形变化，分析频域波形。根据时域、频域的变化验证时域抽样定理 X(ej2肝)=工Xa[j2(fkfs)]o改变Tk程序中频率轴变量为：f(Hz)模拟频率、Q(rad/s)模拟角频率和co(rad)数字频率，观察频率特性的周期变化，并用所学理论对以上各种情况加以解释。例如图1.4给出了单边减幅余弦抽样信号x(n)的时域和幅频特性曲线，由此图可以观察分析选取抽样频率fs=2fmax=2X150=300Hz时的时域和幅频特性。SampleSignalSampleSignalFourierTransform图1.4xa⑴的时域和幅频特性曲线(3)系统响应分析①生成实验用的输入序列x(n)和系统单位冲激响应序列h(n)输入序列：x(n)=Rio(n)单位冲激响应序列:h(n)=沏)+2.5即-1)+2.5<n-2)+即-3)②时域离散信号、系统和系统响应分析观察系统h(n)对输入信号x(n)=Rio(n)的响应特性y(n)。即利用线性卷积求系统响应y(n),调用有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。可以直接调用MATLAB语言中的卷积函数convoconv用于两个有限长度序列的卷积，它假定两个序列都从 n=0开始。调用格式如下：y=conv(x,h)改变输入信号x(n)长度，打印输出y(n)图形并判断y(n)图形及其非零值序列长度是否与理论计算结果一致。③卷积定理的验证a.调用序列博里叶变换数值计算子程序，求输出 y(n)的频率响应Y(ej刃，观察|Y(ej3k)|特性曲线。b.调用序列博里叶变换数值计算子程序，分别求输入序列 x(n)的频率响应X(ej力和离散系统的冲激响应序列h(n)的频率响应H(ejWk),计算Y(ej"尸X(ej，xH(ej，，并绘出|Y(d3k)|曲线。c.将a.与b.分别计算的幅频特性曲线|Y(d3k)|进行比较，验证时域卷积定理。图1.5给出系统响应分析实验程序流程框图，供实验者参考。开始产生输入信号序列x(n)和单位冲激响应序列h(n)

计算离散系统输出信号厂列y(n)=x(n)*h(n)

用序列傅氏变换数值计算，求Y(ejak)

分别绘制y(n),|Y(d力|图形

用序列傅氏变换数值计算求X(ej3k)并绘制其图形用序列傅氏变换数值计算，求H(ej3k)并绘制其图形

求Y(ej，=X(ej3k)XH(ej3k)绘制y(n)的Y(e*图形与丫得业)=X(e*XH(ejak)的图形—结束

图1.5系统响应分析实验的程序框图.思考题在分析理想抽样序列特性的实验中，采样频率不同时，相应理想抽样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？.实验报告要求(1)简述实验目的及实验原理。(2)按实验要求附上实验过程中编制的 MATLAB源程序及实验中涉及的所有信号序歹h系统单位冲激响应及系统响应序列的时域和幅频特性曲线， 并对所得结果进行分析和解释。(3)结合所学理论知识总结实验中的主要结论。(4)简要回答思考题。实验二：用FFT做谱分析.实验目的(1)进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质)。(2)熟悉FFT算法原理和FFT子程序的应用。(3)学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。.实验步骤(1)复习DFT的定义、性质和用DFT作谱分析的有关内容。(2)复习按时间抽选法FFT算法原理及相应的运算流图(3)编制信号产生子程序，产生以下典型信号供谱分析用：x1(n)=R4(n)n1,0n3x2(n)=8n,4n70 ,其他n4 n, 0 n 3x3(n)= n 3, 4 n 70 , 其他nx4(n)=cos(m/4)x5(n)=sin(m/8)x6(t)=cos8疝+cos16d+cos20疝应当注意，如果给出的是连续信号 xa(t),则首先要根据其最高频率确定抽样频率fs以及由频率分辨率选择抽样点数 N,然后对其进行软件抽样(即计算x(n尸xa(nT),0<n<N-1),产生对应序列x(n)0对信x6(t),频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。请实验者根据 DFT的隐含周期性思考这个问题。(4)编写主程序。图2.1给出了主程序框图，供参考,开始读入长度N调用信号产生子程序产生实验信号调用绘图子程序(函数)绘制时间序列波形图调用FFT子程序(函数)计算信号的DFT调用绘图子程序(函数)绘制|X(k)|曲线C结束)图2.1主程序框图(5)按实验内容要求，上机实验，并写出实验报告。.上机实验内容对2中所给出的信号逐个进行谱分析。下面给出针对各信号的 FFT变换区间N以及对连续信号X6(t)的抽样频率fs,供实验时参考。xi(n),X2(n),X3(n),X4(n),X5(n):N=8,16X6(t)：fs=64(Hz),N=16,32,64.思考题(1)在N=8时，X2(n)和X3(n)的幅频特性会相同吗？为什么？N=16呢？(2)如果周期信号的周期预先不知道，如何用 FFT进行谱分析？.实验报告要求(1)简述实验原理及目的。(2)结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析 说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。(3)总结实验所得主要结论。(4)简要回答思考题。实验三：用双线性变换法设计IIR数字滤波器.实验目的(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。(2)掌握数字滤波器的计算机仿真方法。(3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。.实验内容(1)用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2九时，最大衰减小于1dB;在阻带内［0.34nt］频率区间上，最小衰减大于15dB。(2)以0.02冗为抽样间隔，打印出数字滤波器在频率区间［0,N2］上的幅频响应特性曲线。(3)用所设计的滤波器对实际心电图信号抽样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。.实验步骤(1)复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计 IIR数字滤波器的内容，按照“2.实验内容(1)”的要求设计满足指标的数字滤波器函数 H(z)0现给出满足本实验要求的数字滤波器系统函数:H(z)=(11.2686zH(z)=(11.2686z0.0007378(1z1)62 Z1 ""~~2 Z 1Z 2-0.7051z)(11.0106z 0.3583z)(10.9044z 0.2155z)3=3=Hk(z)k1(3.1)式中Hk(z)=式中Hk(z)=A(12z1z2)Z1Z21Bkz1Ckz2(3.2)A=0.09036B1=1.2686 , C1=—0.7051B2=1.0106 , C2=—0.3583B3=0.9044 , C3=—0.2155由(3.1)式和(3.2)式可见，滤波器H(z)由三个二阶滤波器H1(z)、H2(z)和H3(z)级联组成，如图3.1所示。图3.1滤波器H(z)的组成(2)编写滤波器仿真程序，计算H(z)对心电图信号抽样序列x(n)的响应序列y(n)。设yk(n)为第k级二阶滤波器Hk(z)的输出序列，yk-1(n)为输入序列，如图3.1所示。由(3.2)式可得到差分方程：yk(n)=Ayk-1(n)+2Ayk-1(n-1)+Ayk-1(n-2)+Bkyk(n-1)+Ckyk(n-2) (3.3)当k=1时，yk-1(n)=x(n)。所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代算法得到。即依次对k=1,2,3,求解差分方程(3.3),最后得到y3(n)=y(n)。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用 MATLABfilter函数实现仿真。(3)在通用计算机上运行仿真滤波程序，并调用通用绘图子程序，完成“ 2.实验内容(2)和(3)本实验要用的MATLAB绘图函数见附录。.思考题用双线性变换法设计数字滤波器过程中，变换公式中T的取值，对设计结果有无影响？为什么？.实验报告要求(1)简述实验目的及原理。(2)由所打印的|H(e号|特性曲线及设计过程简述双线性变换法的特点。(3)对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用(4)简要回答思考题。.心电图信号抽样序列x(n)

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号抽样序列样本实验四：用窗函数法设计FIR数字滤波器实验四：用窗函数法设计FIR数字滤波器{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,812,12,10，6,6,6,4,0，0，00,0,-2,-4,0,0，0,-2,-2,00,-2,-2,-2,-2,0"x(n),其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分.实验目的(1)掌握用窗函数法设计FIR数字滤波器的原理和方法。(2)熟悉线性相位FIR数字滤波器特性。了解各种窗函数对滤波特性的影响。.实验原理与方法如果所希望的滤波器的理想频率响应函数为 Hd(ej),则其对应的单位脉冲响应为1 jjhd(n)=—Hd(e)ed (4.1)窗函数设计法的基本原理是用有限长单位脉冲响应序列 h(n)逼hd(n)0由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数。 w(n)将hd(n)截断，并进行加权处理，得到：h(n)=hd(n)w(n) (4.2)h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数 H(ej3)为N1(4.3)H(ej)=h(n)e(4.3)n0式中，N为所选窗函数w(n)的长度。由第七章可知，用窗函数法设计的滤波器性能取决于窗函数 w(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度No各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见第这样选定窗函数类型和长度N后，求出单位脉冲响应h(n)=hd(n)-w(n),并按式(4.3)求出H(ej3)。H(ej今是否满足要求，要进行验算。一般在h(n)尾部加零使长度满足2的整数次幕，以使用FFT计算H(ejw)。如果要观察细节，补零点数增多即可。如果H(ejw)不满足要求，则要重新选择窗函数类型和长度N,再次验算，直至满足要求。如果要求线性相位特性，则h(n)还必须满足：h(n)=±h(N-1-n)根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性，可选择h(n)=h(N—1—n)一类，而不能选h(n)=-h(N-1-n)一类。3.实验内容及步骤

(1)复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。(2)编写程序。①编写能产生矩型窗、升余弦窗、改进开余弦窗和二阶开余弦窗的窗函数子程序 。②编写主程序。主程序框图如图4.1所示，仅供参考。其中幅度特性要求用dB表示。读入窗口长度N读入窗口长度N计算hd(n)调用窗函数子程序求w(n)计算h(n)=hd(n)w(n) * 调用子程序(函数)计算H(k)=DFT[h(n)]调用绘图子程序(函数)绘制|H(k)|幅度相位曲线C结束)

图4.1用窗函数法设计滤波器主程序框图H(k)=DFT[h(n)]H(k)=HR(k)+jHi(k)|H(k)|=jHR(k)H：(k)2回图时，用20lg|H(k)|打印幅度特性。第k点对应的频率ok=2-ko为使曲线包络更接近H(d)的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h(n)末尾补零的方法，使长度变为64,再进行64点DFT,则可得到更精确的幅度衰减特性曲线。FFT子程序可调用MALTAB函数fft实现。(3)上机实验内容。①用开余弦窗设计一线性相位低通FIR数字滤波器，截止频率a=—rad。窗口长度4N=15,33。要求在两种窗口长度情况下，分别求出h(n),打印出相应的幅频特性和相频特性曲线，观察3dB带宽20dB带宽。总结窗口长度N对滤波特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数 Hd(d)即HdHd(ej)0,其中 =(N-1)/2jjn 1jjn sincnhd(n)——Hd(e)ed=——eed= 2 nc②n=33,心=N4,用四种窗函数设计线性相位低通滤波器。 绘制相应的幅频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。.思考题(1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减， 如何用窗函数法设计线性相位低通滤波器？写出设计步骤。(2)如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为 ⑴1和试求理想带通的单位脉冲响应hd(n)。.实验报告要求(1)简述实验目的及原理。(2)按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波特性的影响。(3)总结用窗函数法设计FIR滤波器的主要特点。(4)简要回答思考题。附录：MATLAB简介MATLAB(ANTrixLABoratory)是Mathwords公司于1982年推出的一套高性能的数值计算和可视化软件。它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个方便且界面友好的用户环境。 MATLAB的推出得到了各个领域专家学者的广泛关注。其强大的扩展功能为各个领域的应用提供了基础。由各个领域的专家学者相继推出了十几种MATLAB工具箱，信号处理工具箱就是其中之一，借助这一工具箱进行信号处理方面的研究，可大大节省时间。这些工具箱函数的使用可通过Help命令得到。为了便于使用，下面简要介绍MATLAB的启动、编程和运行程序的过程，最后给出实验中可能用到的工具箱函数。(考虑到上机方便，附录中外文均采用正白体。)1、启动MATLAB系统及编程在Windows桌面上双击MATLAB快捷方式图标，启动MATLAB系统，Windows工作平台上自动弹出一个MALAB窗口，其中有CommandWindow窗口，称其为MATLAB命令窗口。该窗口是用户与MATLAB解释器进行通信的工作环境，所有的MATLAB函数、命令及程序都