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

《数字信号处理》上机实验指导书 -PAGE1-《数字信号处理》上机实验指导书一、引言“数字信号处理”是一门理论和实验密切结合的课程，为了深入地掌握课程内容，应当在学习理论的同时，做习题和上机实验。上机实验不仅可以帮助学生深入地理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。所以，根据本课程的重点要求编写了四个实验。第一章、二章是全书的基础内容，抽样定理、时域离散系统的时域和频域分析以及系统对输入信号的响应是重要的基本内容。由于第一、二章大部分内容已经在前期《信号与系统》课程中学习完，所以可通过实验一帮助学生温习以上重要内容，加深学生对“数字信号处理是通过对输入信号的一种运算达到处理目的”这一重要概念的理解。这样便可以使学生从《信号与系统》课程顺利的过渡到本课程的学习上来。第三章、四章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(t)进行理想抽样的过程可用（1.1）式表示。=δT(t)——（1.1）其中为xa(t)的理想抽样，δT(t)为周期冲激脉冲，即——（1.2）的傅里叶变换为=——（1.3）（1.3）式表明为的周期延拓，其延拓周期为抽样角频率(Ωs=2π/T)。抽样前后信号的频谱示意图见“参考教材数字信号处理教程（程佩青）图1-29”。只有满足抽样定理时，才不会发生频率混叠失真。数字信号处理教程（程佩青）在计算机上用高级语言编程直接按（1.3）式计算理想抽样的频谱很不方便。下面导出用序列的傅里叶变换来计算的公式。将（1.2）式代入（1.1）式并进行傅里叶变换，==

==——（1.4）式中的xa(nT)就是采样后得到的序列x(n)，即x(n)=xa(nT)x(n)的序列傅里叶变换为X(ejω)=——（1.5）比较（1.5）和（1.4）可知=X(ejω)|ω=ΩT——（1.6）这说明两者之间只在频率度量上差一个常数因子T。实验过程中应注意这一差别。离散信号和系统在时域均可用序列来表示。序列图形给人以形象直观的印象，它可加深我们对信号和系统的时域特征的理解。本实验还将观察分析几种信号及系统的时域特性。为了在数字计算机上观察分析各种序列的频域特性，通常对X(ejω)在[0，2π]上进行M点采样来观察分析。对长度为N的有限长序列x(n)，有——（1.7）其中,k=0,1,…,M-1通常M应取得大一些，以便观察谱的细节变化。取模||可绘出幅频待性曲线。一个时域离散线性非移变系统的输入/输出关系为y(n)=x(n)*h(n)=——（1.8）这里，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.9）(1.9)式右边的相乘是在各频点{ωk}上的频谱值相乘。3.实验内容(1)连续信号分析①连续时间非周期信号的选择1单边指数脉冲；2双边指数脉冲；3钟形脉冲；4抽样脉冲；1.xa(t)=E*exp(-at)(a>0);2.xa(t)=E*exp(-a*abs(t))(a>0)

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

②用MATLAB编制程序开始结束对xa(t)抽样产生序列x(n)求序列x(n)的傅氏变换X(e开始结束对xa(t)抽样产生序列x(n)求序列x(n)的傅氏变换X(ejω)分别绘制x(n)，|X(ej2πfT)|图形确定抽样频率fs,产生离散时间变量t=nT图1.3离散信号分析程序框图③时域观察，频域分析调整抽样频率fs=2fmax、fs>2fmax和fs<2fmax，观察时域波形变化，分析频域波形。根据时域、频域的变化验证时域抽样定理X(ej2πfT)=。改变程序中频率轴变量为：f(Hz)模拟频率、Ω(rad/s)模拟角频率和ω(rad)数字频率，观察频率特性的周期变化，并用所学理论对以上各种情况加以解释。例如图1.4给出了单边减幅余弦抽样信号x(n)的时域和幅频特性曲线，由此图可以观察分析选取抽样频率fs=2fmax=2×150=300Hz时的时域和幅频特性。图1.4xa(t)的时域和幅频特性曲线

(3)系统响应分析①生成实验用的输入序列x(n)和系统单位冲激响应序列h(n)

输入序列：x(n)=R10(n)

单位冲激响应序列：h(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)②时域离散信号、系统和系统响应分析

观察系统h(n)对输入信号x(n)=R10(n)的响应特性y(n)。即利用线性卷积求系统响应y(n)，调用有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。可以直接调用MATLAB语言中的卷积函数conv。conv用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。

调用格式如下：y=conv(x,h)

改变输入信号x(n)长度，打印输出y(n)图形并判断y(n)图形及其非零值序列长度是否与理论计算结果一致。③卷积定理的验证

a.调用序列博里叶变换数值计算子程序，求输出y(n)的频率响应Y(ejωk),观察|Y(ejωk)|特性曲线。

b.调用序列博里叶变换数值计算子程序，分别求输入序列x(n)的频率响应X(ejωk)和离散系统的冲激响应序列h(n)的频率响应H(ejωk)，计算Y(ejωk)=X(ejωk)×H(ejωk)，并绘出|Y(ejωk)|曲线。

c.将a.与b.分别计算的幅频特性曲线|Y(ejωk)|进行比较，验证时域卷积定理。

图1.5给出系统响应分析实验程序流程框图，供实验者参考。开始结束产生输入信号序列x(n)和单位冲激响应序列h(n)用序列傅氏变换数值计算，求X(ejω开始结束产生输入信号序列x(n)和单位冲激响应序列h(n)用序列傅氏变换数值计算，求X(ejωk)并绘制其图形分别绘制y(n)，|Y(ejωk)|图形用序列傅氏变换数值计算，求H(ejωk)并绘制其图形求Y(ejωk)=X(ejωk)×H(ejωk)计算离散系统输出信号序列y(n)=x(n)*h(n)用序列傅氏变换数值计算，求Y(ejωk)绘制y(n)的Y(ejωk)图形与Y(ejωk)=X(ejωk)×H(ejωk)的图形图1.5系统响应分析实验的程序框图4.思考题在分析理想抽样序列特性的实验中，采样频率不同时，相应理想抽样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？5.实验报告要求(1)简述实验目的及实验原理。(2)按实验要求附上实验过程中编制的MATLAB源程序及实验中涉及的所有信号序列、系统单位冲激响应及系统响应序列的时域和幅频特性曲线，并对所得结果进行分析和解释。(3)结合所学理论知识总结实验中的主要结论。(4)简要回答思考题。实验二：用FFT做谱分析1.实验目的(1)进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。(2)熟悉FFT算法原理和FFT子程序的应用。(3)学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。2.实验步骤(1)复习DFT的定义、性质和用DFT作谱分析的有关内容。(2)复习按时间抽选法FFT算法原理及相应的运算流图(3)编制信号产生子程序，产生以下典型信号供谱分析用：x1(n)=R4(n)x1=ones(1,4)=[1,1,1,1]x2(n)=x2=[1,2,3,4,4,3,2,1…0]x3(n)=x4(n)=cos(π/4)nx5(n)=sin(π/8)nx6(t)=cos8πt+cos16πt+cos20πt应当注意，如果给出的是连续信号xa(t)，则首先要根据其最高频率确定抽样频率fs以及由频率分辨率选择抽样点数N，然后对其进行软件抽样(即计算x(n)=xa(nT)，

0≤n≤N-1)，产生对应序列x(n)。对信x6(t)，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。请实验者根据DFT的隐含周期性思考这个问题。

开始结束开始结束读入长度N调用信号产生子程序产生实验信号调用绘图子程序(函数)绘制时间序列波形图调用FFT子程序(函数)计算信号的DFT调用绘图子程序(函数)绘制|X(k)|曲线

图2.1主程序框图(5)按实验内容要求，上机实验，并写出实验报告。3．上机实验内容(1)对2中所给出的信号逐个进行谱分析。下面给出针对各信号的FFT变换区间N以及对连续信号x6(t)的抽样频率fs，供实验时参考。x1(n),x2(n),x3(n),x4(n),x5(n)：N=8,16x6(t)：fs=64(Hz),N=16,32,64(2)令x(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换，X(k)=DFT[x(n)]并根据DFT的对称性，由X(k)求出X4(k)=DFT[x4(n)]和X5(k)=DFT[x5(n)]，并与(1)中所得结果比较。［提示：取N=16时，x4(n)=x4(N-n),x5(n)=-x5(N-n)。］(3)令x(n)=x4(n)+jx5(n)，重复(2)。4.思考题(1)在N=8时，x2(n)和x3(n)的幅频特性会相同吗？为什么？N=16呢？(2)如果周期信号的周期预先不知道，如何用FFT进行谱分析？5.实验报告要求(1)简述实验原理及目的。(2)结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。(3)总结实验所得主要结论。(4)简要回答思考题。实验三：用双线性变换法设计IIR数字滤波器1.实验目的(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。(2)掌握数字滤波器的计算机仿真方法。(3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。2.实验内容(1)用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π,π]频率区间上，最小衰减大于15dB。(2)以0.02π为抽样间隔，打印出数字滤波器在频率区间[0,π/2]上的幅频响应特性曲线。(3)用所设计的滤波器对实际心电图信号抽样序列（在本实验后面给出）进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。3.实验步骤(1)复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，按照“2.实验内容(1)”的要求设计满足指标的数字滤波器函数H(z)。现给出满足本实验要求的数字滤波器系统函数：H(z)=

=——（3.1）

式中Hk(z)=,k=1,2,3——（3.2）A=0.09036

B1=1.2686,C1=－0.7051

B2=1.0106,C2=－0.3583

B3=0.9044,C3=－0.2155由(3.1)式和(3.2)式可见，滤波器H(z)由三个二阶滤波器H1(z)、H2(z)和H3(z)级联组成，如图3.1所示。H1H1(z)H2(z)H3(z)x(n)y1(n)y2(n)y3(n)=y(n)H(z)图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绘图函数见附录。4.思考题用双线性变换法设计数字滤波器过程中，变换公式s=中T的取值，对设计结果有无影响？为什么？5.实验报告要求(1)简述实验目的及原理。(2)由所打印的|H(ejω)|特性曲线及设计过程简述双线性变换法的特点。(3)对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。(4)简要回答思考题。6.心电图信号抽样序列x(n)人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号抽样序列样本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，｝实验四：用窗函数法设计FIR数字滤波器1.实验目的(1)掌握用窗函数法设计FIR数字滤波器的原理和方法。(2)熟悉线性相位FIR数字滤波器特性。(3)了解各种窗函数对滤波特性的影响。2.实验原理与方法如果所希望的滤波器的理想频率响应函数为Hd(ejω)，则其对应的单位脉冲响应为hd(n)=——（4.1）窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数。w(n)将hd(n)截断，并进行加权处理，得到：h(n)=hd(n)w(n)——（4.2）h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数H(ejω)为H(ejω)=——（4.3）式中，N为所选窗函数w(n)的长度。由第七章可知，用窗函数法设计的滤波器性能取决于窗函数w(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见第七章。这样选定窗函数类型和长度N后，求出单位脉冲响应h(n)=hd(n)·w(n)，并按式(4.3)求出H(ejω)。H(ejω)是否满足要求，要进行验算。一般在h(n)尾部加零使长度满足2的整数次幂，以便用FFT计算H(ejω)。如果要观察细节，补零点数增多即可。如果H(ejω)不满足要求，则要重新选择窗函数类型和长度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调用子程序(函数)计算H(k)=DFT[h(n)]结束计算hd(n)调用绘图子程序(函数)绘制|H(k)|幅度相位曲线调用窗函数子程序求w(n)计算h(n)=hd(n)w(n)设H(k)=DFT[h(n)]H(k)=HR(k)+jHI(k)|H(k)|=画图时，用20lg|H(k)|打印幅度特性。第k点对应的频率ωk=。为使曲线包络更接近H(ejω)的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h(n)末尾补零的方法，使长度变为64，再进行64点DFT，则可得到更精确的幅度衰减特性曲线。FFT子程序可调用MALTAB函数fft实现。(3)上机实验内容。①用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率ωc=rad。窗口长度N=15,33。要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频特性曲线，观察3dB带宽20dB带宽。总结窗口长度N对滤波特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数Hd(ejω)，即其中=(N-1)/2==②n=33，ωc=π/4，用四种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。4.思考题(1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出设计步骤。(2)如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd(n)。5.实验报告要求(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函数、命令及程序都要在该命令窗口下运行。在命令窗口中，用户可以发出MATLAB命令。每条命令输入并接回车键后，MATLAB系统便解释并执行之，显示命令执行结果。如果不想显示中间结果，可在相应命令后加分号即可。应当注意，在命令窗口输人命令只适合一些简单的运算程序，不便于编辑修改较复杂的程序。这是因为直接在命令窗口输入的命令一经执行后，只记录其运行结果变量，如果有错误，就得再输入一遍（当然，有较简便的重新输入方法）。这样一来调试程序太麻烦。M文件为这种编辑运行程序过程提供了方便。在MATLAB系统中，根据功能可将MATLAB系统所使用的外部文件分成几类，并用不同的扩展名作为标识。(1)M文件以字母m为扩展名，一般为ASCII码文本文件，可以用任何文本编辑器进行编辑。(2)MAT文件是MATLAB系统的二进制数据文件，用于保存MATLAB系统所使用的数据。(3)MEX文件是经过MATLAB编译系统编译的函数二进制可执行文件，可直接调入MATLAB系统中运行，运行速度比M文件快得多，所以常把已经调试好且比较大的M文件编译成MEX文件，供以后使用。M文件又分为程序M文件和函数M文件。下面简要介绍实验过程中要用到的程序M文件的编辑与运行方法。如前述，M文件为ASCII文本文件，可以用任何文本编辑器进行编辑。本节只介绍在MATLAB命令窗口状态下M文件的编辑运行与修改。1)新建M文件在MATLAB命令窗口状态，选File-New-Mfile，MATLAB系统弹出文件编辑窗口，类似于Windows的记事本（低版本MATLAB就是调用Windows记事本），这时一切编辑就如同记事本一样。2)运行与修改M文件MATLAB采用的是行命令模式，用户每输入一条命令，回车后MATLAB就解释并执行这条命令，再根据要求显示运算结果。此外，MATLAB也可执行某个文件中的MATLAB语句序列（类似于DOS的批处理文件）。这两种模式一起构成了MATLAB的解释环境。当一个M文件编辑存盘后，在MATLAB命令窗口输入M文件名并回车，则MATLAB系统就逐行解释并执行该M文件中的命令序列（即程序）。若发现错误，则打开M文件修改，再运行。应当注意，如果已编辑好的M文件不在MATLAB系统搜索路径中，用户应当用path命令在MATLAB的搜索路径中添加新的搜索路径。例如，M文件存在C：\MYFILES路径，下面的命令：path(‘C:MYFILES’，path)；将搜索路径改为在搜索完当前目录之后，先搜索C:\MYFILES，再在当前的搜索路径目录中搜索。2、实验用MATLAB函数简介为了便于实验时查阅，下面简要介绍几种实验中可能要用到的MATLAB函数。1）abs功能：求绝对值（模值）。格式：y=abs(x)说明：y=abs(x)用于计算x的绝对值，当x为复数时，得到的是复数的模值。当x为字符串时，abs(x)得到字符串的各个字符的ASCII码，例如，x=’123’，则abs(x)得到：4950512）angle功能：求相角。格式：φ=angle(h)说明：φ=angle(h)用于求复矢量或复矩阵的相角（以弧度为单位），相角介于-π和+π之间。例如，对复数h可用两种方法表示：h=x+jy=mejφ则m和φ可由x+jy表示格式的复数求得m=abs(h)φ=angle(h)当然，也可由m和φ求取h=x+jy格式中的x和y：h=m.*exp(i*φ)x=real(h)y=imag(h)3）conv功能：求卷积。格式：c=conv(a,b)说明：conv(a,b)用于求矢量a和b的卷积，即c(n)=a(k+1)b(n-k),n=1,2,…式中N为矢量a和b的最大长度。例如，当a=[123],b=[156]时，则c=conv(a,b)c=413282718此函数可直接用于求两个有限长序列的卷积。设x(n)和h(n)的长度分别为M和N，则计算二者卷积的MATLAB语句如下：y=conv(x,h)y的长度为N+M-1。4)filter功能：利用IIR滤波器或FIR滤波器对数据进行滤波。格式：y=filter(b,a,x)[y,zf]=filter(b,a,x)y=filter(b,a,x,zi)说明：filter利用数字滤波器对数据进行滤波，其实现采用直接五型结构，因而适用于IIR和FIR两种滤波器。滤波器的系统函数为即滤波器系数a=[a0a1a2…aN],b=[b0b1…bM]，输入序列矢量为X。这里的标准形式为a0=1，如果输入矢量a时，a0≠1，则MATLAB将自动进行归一化系数的操作；如果a0=0，则给出出错信息。y=filter(b,a,x)利用给定系数矢量a和b对x中的数据进行滤波，结果放入y矢量中，y的长度取max(N,M)。y=filter(b,a,x,zi)可在zi中指定x的初始状态。[y,zf]=filter(b,a,x,zi)除了得到矢量y外，还得到x的最终状态矢量zf。例C1设时域离散系统可由如下差分方程描述：y(n)=y(n-1)-0.9y(n-2)+x(n),n(1)计算n=0,…，120上的单位脉冲响应；(2)计算n=0,…，120上的单位阶跃响应。解可利用filter函数直接求解。MATLAB程序为exc1.m。%Examplecla=[1,-1,0.9];b=1;x=zeros(1,120);n=[0:119];x(1)=1;%x(n)=δ(n)h=filter(b,a,x);subplot(2,1,1);stem(n,h,'.');title('ImpulseResponse');xlabel('n');ylabel('h(n)');x=ones(1,120);%x(n)=u(n)s=filter(b,a,x);subplot(2,1,2);stem(n,s,'.');title('StepResponse');xlabel('n');ylabel('s(n)');系统的单位脉冲响应h(n)和单位阶跃响应s(n)如图C．1所示。仿照该例题，可调用filter函数进行实验三中的IIR滤波仿真。

图C．1系统的单位脉冲响应h(n)和阶跃响应s(n)5)freqz功能：数字滤波器的频率响应。格式：[h,w]=freqz(b,a,n)[h,f]=freqz(b,a,n,Fs)h=freqz(b,aw)h=freqz(b,a,f,Fs)freqz(b,a)说明：freqz用于计算数字滤波器的频率响应函数H(ejω)。H(z)具有如下形式:H(z)=矢量a=[1a1a2…aN]和b=[b0b1…bM]分别表示分母多项式A(z)和分子多项式B(z)的系数。[h,w]=freqz(b,a,n)可得到数字滤波器的n点频率响应值，这n个点均匀地分布在[0，π]上，并将这n个频点的频率记录在w中，相应的频响值记录在h中。要求n为大于零的整数，最好为2的整数次幂，以便采用FFT计算，以提高速度。[h,f]=freqz(b,a,n,Fs)用于对H(ejω)在［0，Fs/2］上等间隔采样n点，采样点频率及相应频响值分别记录在f和h中。由用户指定Fs（以Hz为单位）值。h=freqz(b，a，w)用于对H(ejω)在[0，2π]上进行采样，采样频率点由矢量w指定。h=freqz(b,a,f,Fs)用于对H(ejω)在[0，Fs]上采样，采样频率点由矢量f给定。freqz(b,a,n)用于在当前图形窗口中绘制出幅频和相频特性曲线。6)impz功能：计算H(z)相应的单位脉冲响应h(n)。格式：[h,t]=impz(b,a)[h,t]=impz(b,a,z)[h,t]=impz(b,a,n,Fs)impz(b,a)其中a和b分别为H(z)的分母和分子多项式系数。该函数的详细说明请看MATLABHelp命令。7）fft功能：一维快速傅里叶变换(FFT)。格式：y=fft(x)y=fft(x,n)说明：fft函数用于计算矢量或矩阵的离散傅里叶变换。y=fft(x)利用FFT算法计算矢量x的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。当x长度为2的整数次幂时，fft采用基2FFT算法，否则采用稍慢的混合基算法。y=fft(x,n)采用n点FFT。当x长度小于n时，fft函数自动在x尾部补零，以构成n点数据；当x的长度大于n时，fft截取x的前面n点数据进行FFT。8）ifft功能：一维逆快速傅里叶变换（IFFT）。格式：y=ifft(x)y=ifft(x,n)9)plot功能：线型绘图函数。格式：plot(v)plot(x,y)说明：plot(v)是一种最简单的调用方式，v是长度为n的数值向量。Plot(v)的作用是在坐标系中顺序地用直线连接顶点{i,v(i),i=1,2,…,n}生成一条折（曲）线。当向量元素充分多时，既可生成一条光滑的曲线。在实验中，若FFT点数足够多时，用plot打印的幅频特性就很接近x(ejω)连续曲线。Plot(x,y)称为参数式，参数x和y都是长度为n的向量。plot(x,y)在坐标系中生成顺序连接顶点{x(i),y(i),i=1,2,…,n}的折（曲）线，这种调用可被用来生成参数方程的图形。10）stem功能：绘制离散序列图。格式：stem(y)stem(x,y)stem(…,’线端符号’)stem(…,’线型’)stem(…,’线型’,’线端符号’)说明：stem(y)和stem(x,y)分别与plot(x)和plot(x,y)的绘图规则相同，只是stem绘制的是离散序列图（或称为“杆状”图），序列线端为圆圈。实验中用于绘制时域序列x(n)的波形图和序列的离散傅里叶变换x(k)的幅度图。后面三种调用方式中的‘线型’用于指定图中序列线的线型和颜色，“线端符号”用于指定序列线端点表示序列值的绘图形状。线型、颜色和线端符号如表C.1所示。表C.1线型、颜色与线端符号线型颜色线端符号类型符号类型符号类型符号实线(缺省)-红r实点.点线:绿g星号*点虚线-.蓝b圆圈o虚线--黑k三角形^表中只列出几种共实验用，其他更多的符号和颜色清查阅MATLAB书。11）subplot功能：多坐标设置与定位当前坐标系。格式：subplot(m,n,k)说明：s