《数字信号处理》2

数字信号处理实 验 指 导 书周春 编 写适用专业： 电子信息工程 通信工程 安徽建筑工业学院电子与信息工程学院2008年 09月前 言数字信号处理是一门理论和实践密切结合的课程，为深入掌握课程内容，应当在学习理论的同时上机实验。上机实验不仅可以帮助学生深入的理解和消化基本理论，而且能锻炼初学者独立解决问题的能力。本书根据课程重点编写了四个实验。采样定理，时域离散系统的时域和频域分析以及系统对输入信号的响应是数字信号处理的基本内容，在学习完该部分内容后，进行实验一，以帮助读者掌握以上内容。同时加深学生对“数字信号处理是通过对输入信号的一种运算达到处理目的的一种运算”的理解DFT、FFT是数字信号处理的重要数学工具，它具有广泛的使用内容。学习完该部分内容后，进行实验二“用FFT对信号进行谱分析”，通过该实验加深理解DFT的基本概念、基本性质，及其快速算法FFT。DFT和FFT是学习数字信号处理必须掌握的内容。数字滤波器的基本理论和设计方法是数字信号处理技术的重要内容。学习这一部分时，应重点掌握IIR和FIR两种滤波器的基本设计方法。IIR滤波器的单位脉冲响应是无限长的，设计方法是先设计模拟滤波器，然后通过S-Z平面转换，求出响应滤波器的系统函数。这里的平面转换有两种方法，即脉冲响应不变法和双线性变换法。后者没有频率混叠的缺点，且转换简单，是一种普遍应用的方法。FIR滤波器的单位脉冲响应是有限长的，设计滤波器的目的即是求出符合要求的单位脉冲响应。窗函数法是一种基本的，也是重要的设计方法。数字信号处理实验适用于电子信息工程、通信工程等专业学生的实践性教学环节，通过该课程的学习，使学生能够熟练掌握数字信号处理的一般方法，熟练使用Matlab作为编程工具，实现在课堂上学习过的各种数字信号处理方法与滤波器设计方法，在此基础上养成严谨认真的科学工作作风，为今后的工作打下良好的基础。在本课程的学习中，学生应达到如下要求：1、提前预习课本、实验指导书和参考书籍，巩固和加深对实验中要用到基本理论知识的理解；2、在实验过程中认真观察分析实验结果,运用所学知识解决实验中发生的问题； 3、能熟练使用编程工具Matlab； 4、正确处理实验数据，并完成符合要求的实验报告。本实验指导书共包括四个实验，分别是：信号、系统及系统响应、用FFT作谱分析、用双线性变换法设计IIR数字滤波器、用窗函数法设计FIR数字滤波器，均为必做实验。主讲教师可以根据实际教授的专业和课时情况适当增、减部分实验内容。学生在做完每个实验后，必须独立写出详细的实验报告，包括实验方法、实验过程（步骤）及结果、实验结果分析等等。实验简要介绍序号实验项目名称内 容 提 要主要仪器设备及台套数实验类型实验要求计划学时每组人数1信号、系统及系统响应1、用matlab画出指定的连续信号，并对其采样，观察不同采样频率的采样结果。2、用matlab观察信号的时域、频域响应，验证卷积定理计算机Matlab验证必做2402用FFT作谱分析1、用matlab显示指定信号并做谱分析2、观察不同点数FFT结果计算机Matlab验证必做2403用双线性变换法设计IIR数字滤波器用双线性变换法设计巴特沃斯低通IIR数字滤波器，显示其幅频率响应，并用其对心电信号滤波计算机Matlab验证必做2404用窗函数法设计FIR数字滤波器编写程序实现矩形窗、升余弦窗、改进升余弦窗和二阶升余弦窗的FIR数字滤波器计算机Matlab验证必做240合计8目 录前 言1实验一：信号、系统及系统响应5实验二：用FFT作谱分析11实验三：用双线性变换法设计IIR数字滤波器16实验四：用窗函数法设计FIR数字滤波器19实验一：信号、系统及系统响应 实验学时：2实验类型：验证 实验要求：必修一、实验目的 1、熟悉连续信号经理想采样前后的频谱变化关系，加深对采样定理的理解。2、熟悉时域离散系统的时域特性。3、利用卷积方法观察分析系统的时域特性。4、掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。二、实验内容1、认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容。2、编制实验用主程序及相应子程序，包括（1）信号产生子程序，用于产生实验中要用到的下列信号序列：a、采样信号序列：对下面的连续信号： （1.1）进行采样，可得到采样序列 （1.2）其中A为幅度因子，a为衰减因子，0是模拟角频率，T为采样间隔。这些参数都要在实验过程中由键盘输入，产生不同的xa(t)和xa(n)。b、单位脉冲序列： c、矩形序列： （2）系统单位脉冲响应序列产生子程序。本实验要用到两种FIR系统a、b、（3）有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。可以直接调用MATLAB语言中卷积函数conv。conv用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。调用格式如下：y = conv(x,h)其中参数x和y是两个已赋值的行向量序列。在完成编制上述子程序的基础上，编制本实验主程序。图1.1给出主程序流程框图，供参考。图1.1三、实验原理、方法和手段采样是连续信号数字处理的第一个关键环节。对采样过程的研究不仅可以了解采样前后信号时域和频域发生的变化以及信号信息不丢失的条件，而且可以加深对傅里叶变换、Z变换和序列傅里叶变换之间关系式的理解。我们知道，对一个连续信号xa(t)进行理想采样的过程可以用下式表示 （1.3）其中为的理想采样，p(t)为周期冲激脉冲，即 （1.4）的傅里叶变换为 （1.5）上式表明，为的周期延拓，其延拓周期为采样角频率（s=2/T）。只有满足采样定理时，才不会发生频率混叠失真。在计算机上用高级语言编程，直接按照1.5式计算的频谱很不方便，下面导出用序列的傅里叶变换来计算的公式。将1.4式代入1.3式，并进行傅里叶变换， （1.6）式中的xa(nT)就是采样后的序列x(n)，即：x(n)= xa(nT)x(n)的傅里叶变换为 （1.7）比较1.7和1.6可知 （1.8）说明两者之间只在频率度量上差一个常数因子T。实验过程中应注意这一差别离散信号和系统在时域均可以用序列来表示，序列图形给人以形象直观的印象，它可以加深我们对信号和系统的时域特征的理解。本实验还将观察分析几种信号及系统的时域特性。为了在数字计算机上观察分析各种序列的频域特性，通常对X(ej)在0，2上进行M点采样来观察分析。对长度为N的有限长序列x(n)，有 （1.9）其中通常M应取的大一些，以便观察谱的细节变化。取模可绘出幅频特性曲线。一个时域离散线性非移变系统的输入/输出关系为 （1.10）这里y(n)为系统的输出序列，x(n)为输出序列。h(n)、x(n)可以是无限长，也可以是有限长，为了计算机绘图观察方便，主要讨论有限长情况。如果h(n)和x(n)的长度分别为N和M，则y(n)的长度为L=N+M-1。这样（1.10）式所描述的卷积运算就是序列移位、相乘和累加的过程，所以编程十分简单。上述卷积运算也可以在频域实现 （1.11）1.11式右边的相乘是在个频点k上的频谱值相乘。四、实验组织运行要求根据本实验的特点、要求和具体条件，采用集中授课形式为主的方式进行。五、实验条件1、计算机2、matlab 图1.2 |xa(jf)|六、实验步骤 1、分析采样序列的特性，产生采样信号序列xa(n)，使。图1.2给出了连续信号xa(t)的幅频特性曲线。由此图可以确定对xa(t)应采用的采样频率。 a、取采样频率fs=1kHz，即T=1ms。观察所得采样xa(n)的幅频特性和图1.2中的在折叠频率附近有无明显差别。应当注意，实验中所得频率是用序列xa(n)的傅里叶变换公式求得的，所以在频率度量存在关系：，为数字频率，为模拟频率。b、改变采样频率，fs=300Hz，观察，并做记录（打印曲线）；进一步降低采样频率fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录（打印）这时的曲线。2、时域离散信号、系统和系统响应分析。a、观察信号xb(n)和系统hb(n)的时域和频域特性；利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n)，比较所求响应y(n)和hb(n)的时域和频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。b、观察系统ha(n)对信号xc(n)的响应特性。利用线性卷积求系统响应y(n)，并判断y(n)的图形及其非零值序列长度是否与理论结果一致，对xc(n)= ha(n)=R10(n)，说出一种定性判断y(n)图形正确与否的方法，调用序列傅里叶变换数值计算子程序，求得，观察特性曲线，定性判断结果的正确性。改变xc(n)的长度，取N=5，重复改实验。注意参数变化的影响，说明变化前后的差异，并解释所得结果。3、卷积定理的验证。将实验2中的信号换成xa(n)，使，重复实验2a，打印曲；对主程序做简单修改，按式1.11计算所得结果之间有无差异？为什么？七、思考题1、在分析理想采样序列特性的实验中，采样频率不同时，相应的理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？2、在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得所得结果之间有无差异？为什么？八、实验报告1、简述实验目的及实验原理。2、按照实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时域和幅频特性曲线，并对所得结果进行分析和解释。3、总结实验中的主要结论。4、简要回答思考题。九、其它说明实验二：用FFT作谱分析 实验学时：2实验类型：验证实验要求：必修一、实验目的1、进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）2、熟悉FFT算法原理和FFT子程序的应用3、学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT二、实验内容1、对以下信号逐个进行谱分析 （2.1）下面给出针对各信号的FFT变换区间N以及对连续信号x6(t)的采样频率fs，供实验时参考。x1(n)，x2(n)，x3(n)，x4(n)，x5(n)：N=8，16 x6(t)：fs=64(Hz)，N=16，32，642、令x(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换，X(k)=DFTx(n)并根据DFT的对称性，由X(k)求出X4(k)=DFTx4(n)和X5(k)=DFTx5(n)，并与1中所得的结果进行比较。提示：取N=16时，x4(n) = x4(N-n)，x5(n) = -x5(N-n)3、令x(n)=x4(n)+jx5(n)，重复2。三、实验原理、方法和手段参考教材第三章、第四章。四、实验组织运行要求根据本实验的特点、要求和具体条件，采用集中授课为主的方式进行。五、实验条件1、计算机2、matlab六、实验步骤1、复习DFT的定义、性质和用DFT作谱分析的有关内容。2、复习FFT算法原理与编程思想，并对照DFT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。3、编制信号产生子程序，产生式2.1中的典型信号供谱分析用。应当注意，如果给出的是连续信号xa(t)，则首先要根据其最高频率确定采样速率fs以及由频率分辨率选择采样点数N，然后对其进行软件采样（即计算x(n)= xa(nT)，0nN-1），产生对应序列x(n)。对信号x6 (t)，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。请实验者根据DFT的隐含周期性思考这个问题。4、编写主程序，图2.1给出了主程序框图，供参考。5、按实验内容要求，上机实验，并写出实验报告。本实验的主程序比较简单，直接根据图2.1所给的框图编写程序即可。编程的难点是FFT子程序。不过，各种语言的FFT子程序都可在有关信号处理的程序库中找到。由于C语言当前最普及，所以为了为实验提供方便，下面给出C语言FFT函数，供参考。图2.1/*DITFFT函数（C语言）*/fft基2DITFFT函数要求：指向复数数组指针X，FFT长度为2m，m为正整数FFT输出结果放在输入复数数组中。/*计算N点FFT子程序*/*xr：=信号序列实部，xi：=信号序列虚部，N=FFT变换区间长度N=2M*/*如果信号长度小于N，应该给xr，xi后面补0*/*计算结果X(k)的实部和虚部分别存储在数组xr和xi中*/void Fft(double xr,double xi,int N,int M) int L,B,J,P,k,I; double rPartKB,iPartKB;/*分别代表X（K+B）的实部和虚部*/ double rCf128,iCf128;/*rCf存储旋转因子实部，iCf存储旋转因子虚部*/ /*旋转因子数组长度应根据自己需要调整*/*计算旋转因子*/ double PI2=8.0*aran(1.0) for(i=0;iN,i+) rCfi=cos(i*PI2/N); iCfi=sin(i*PI2/N); ChangeOrder(xr,xi,N);/*调用倒序子程序*/ /*计算各级蝶形*/ for(L=1;L=M;L+) B = (int)(pow(2,(L-1)+0.5); for(J=0;J=B-1;J+) P = J*(int)(pow(2,(M-L)+0.5); rPartKB = xrk+B*rCfP-xik+B*rCfP; iPartKB = xik+B*rCfP+xrk+B*rCfP; xrr+B = xrk rPartKB; xir+B = xik iPartKB; xrk = xrk + rPartKB; xik = xik + iPartKB; /*倒序子程序*/ void ChangeOrder(double xr,double xi,int N) int LH,N1,I,J,K; double T; LH = N/2;J=LH;N1 = N-2; for(I = 1;IN1;I+) if(I=K) J=J-K; K=(int)(K/2+0.5); J=J+K; 七、思考题1、在N=8时，x2(n)和x3(n)的幅频特性会相同吗？为什么？N=16呢？2、如果周期信号的周期预先不知道，如何用FFT进行谱分析呢？八、实验报告1、简述实验原理及目的。2、结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。3、总结实验所得主要结论。4、简要回答思考题。九、其它说明实验三：用双线性变换法设计IIR数字滤波器 实验学时：2实验类型：验证实验要求：必修一、实验目的1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。2、掌握数字滤波器的计算机仿真方法。3、通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。二、实验内容1、用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为在通频带内频率低于0.2时，最大衰减小于1dB；在阻带内0.3, 频率区间上，最小衰减大于15dB。2、以0.02为采样间隔，打印出数字滤波器在频率区间0，/2上的幅频响应特性曲线。3、用所设计的滤波器对实际心电图信号采样序列（在本实验后面给出）进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。三、实验原理、方法和手段参考教材第六章四、实验组织运行要求根据本实验的特点、要求和具体条件，采用集中授课为主的方式进行。五、实验条件1、计算机2、matlab六、实验步骤1、复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器，使其系统函数 (3.1)式中 (3.2) 由3.1和3.2式可见，滤波器H(z)由三个二阶滤波器H1(z)、 H2(z)和 H3(z)级联组成，如下图所示图3.1滤波器H(z)的组成2、编写滤波器仿真程序，计算H(z)对心电图信号采样序列x(n)的响应序列y(n)。设为第k级二阶滤波器Hk(z)的输出序列，为输入序列，如图3.1所示，由3.2式可得到差分方程 (3.3)当k=1时，。所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代算法得到。即依次对k=1，2，3，求解差分方程3.3，最后得到。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用MATLAB filter函数实现仿真3、在通用滤波器上运行仿真滤波程序，利用MATLAB绘图，完成实验内容2、和3七、思考题1、用双线性变换法设计数字滤波器过程中，变换公式中T的数值，对设计结果有无影响？为什么？八、实验报告1、简述实验目的及原理。2、由所打印的特性曲线及设计过程简述双线性变换法的特点。3、对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。4、简要回答思考题。九、其它说明心电图信号采样序列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数字滤波器 实验学时：2实验类型：验证实验要求：必修一、实验目的1、掌握用窗函数法设计FIR数字滤波器的原理和方法。2、熟悉线性相位FIR数字滤波器特性。3、了解各种窗函数对滤波器特性的影响。二、实验内容1、用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率。窗口长度N=15，33。要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N对滤波特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数，即，其中2、，用四种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波特性的影响。三、实验原理、方法和手段如果所希望的滤波器理想频率响应函数为，则其对应的单位脉冲响应 为 (4.1)窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼近hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数(n)将hd(n)截断，并进行加权处理，得到： h(n)= hd(n)(n) (4.2)h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为： (4.3)式中，N为所选窗函数(n)的长度。用窗函数法设计的滤波器性能取决于窗函数(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。这样选定窗函数类型和长度N后，求出单位脉冲响应h(n)= hd(n)(n)，并按式4.3的要求求出。是否满足要求，要进行验算。一般在h(n)尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N，再次进行验算，直至满足要求。如果要求线性相位特性，则h(n)还必须满足：h(n)=h(N-1-n)根据上式中的正、负号和长度N的奇、偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波器特性正确选择其中一类。例如，要设计线性相位低通特性，可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。四、实验组织运行要求根据本实验的特点、要求和具体条件，采用集中授课为主的方式进行。五、实验条件1、计算机2、matlab六、实验步骤1、复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。2、编写程序。(1)编写能产生矩形窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序。(2)编写主程序，主程序框图如图4.1所示，仅供参考，其中幅度特性要求用dB表示。图4.1 用窗函数法设计滤波器主程序框图设画图时，用20lg|H(k)|打印幅度特性。用k点对应的频率。为使曲线包络更接近的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h(n)末尾补零的方法，使长度变为64，再进行64点DFT，则可得到更精确的幅度衰减特性曲线。FFT子程序参考实验一，也可以调用MATLABA函数fft实现。七、思考题1、如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出设计步骤。2、如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为1和2，实求理想带通的单位脉冲响应hd(n)。八、实验报告1、简述实验目的及原理。2、按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波特性的影响。3、总结用窗函数法设计FIR滤波器的主要特点。4、简要回答思考题。九、其它说明22