数字信号处理实验参考书

目 录目录1实验一 信号、系统及系统响应2实验二 应用FFT对信号进行频谱分析6实验三 用双线性变换法设计IIR滤波器11实验四 用窗函数设计FIR滤波器16附录A C语言实现数字信号处理算法21附录B MATLAB的信号表示和处理32附录C MATLAB 下的数字信号处理实现示例56附录D MATLAB 下的数字滤波器设计示例64实验一 信号、系统及系统响应一实验目的1 1 熟悉理想采样的性质，了解信号采样前后的频谱变化，加深对采样定理的理解。2 2 熟悉离散信号何系统的时域特性。3 3 熟悉线性卷积的计算编程方法：利用卷积的方法，观察、分析系统响应的时域特性。4 4 掌握序列傅氏变换的计算机实现方法，利用序列的傅氏变换对离散信号、系统及系统响应进行频域分析。二实验原理 （一）连续时间信号的采样 采样是从连续信号到离散时间信号的过渡桥梁，对采样过程的研究不仅可以了解采样前后信号时域何频域特性发生的变化以及信号内容不丢失的条件，而且有助于加深对拉氏变换、傅氏变换、Z变换和序列傅氏变换之间关系的理解。对一个连续时间信号进行理想采样的过程可以表示为信号与一个周期冲激脉冲的乘积，即： （1-1）其中，是连续信号的理想采样，是周期冲激脉冲（1-2）此时，采样信号的拉氏变换可以表示为：,其中是信号的拉氏变换作为拉氏变换的一种特例，信号理想采样的傅立叶变换为： （1-3）由（1-3）式可知，信号理想采样后的频谱式原来信号频谱的周期延拓，其延拓周期等于采样频率。根据香农定理，如果原信号是带限信号，且采样频率高于原信号最高频率的2倍，则采样后的离散序列不会发生频谱混叠现象。在计算机处理时，不采用时计算信号的频谱，而是利用序列的傅立叶变换计算信号的频谱，定义序列：（1-4）根据Z变换的定义，可以得到序列的Z变换为：（1-5）以代替上式中的Z，就可以得到序列的傅立叶变换：（1-6）式（1-3）和式（1-6）具有如下关系：（1-7）由式可知，在分析一个连续时间信号的频谱时，可以通过取样将有关的计算转化为序列傅立叶变换的计算。（二）有限长序列的分析一般来说，在计算机上不可能、也不必要处理连续的曲线，通常，我们的做饭时只要观察、分析在某些频率点上的值。对于长度为N的有限长序列： （1-8）一般只需要在之间均匀的取M个频率点，计算这些点上的序列傅立叶变换： （1-9）其中，。是一个复函数，它的模就是幅频特性曲线。（三）信号卷积一个线性时不变离散系统的响应可以用它的单位冲激响应和输入信号的卷积来表示： （1-6）根据傅立叶变换和 Z变换的性质，应该有： （1-7）（1-8）上两式告诉我们：可以通过对两个序列的移位、相乘和累加计算信号响应；卷积运算可以在频域通过乘积实现。三实验内容及步骤(一) 编制实验用主程序及相应子程序1信号产生子程序，包括：（1）理想采样信号序列：对信号进行理想采样，可以得到一个理想的采样信号序列，其中为幅度因子，是衰减因子，是频率，为采样周期。根据实验内容的需要，这些参量请设计为在程序运行过程中输入。（2）单位脉冲序列（3）矩形序列, 其中取：N=102系统单位冲激响应序列产生子程序，本实验中用到两种FIR系统：（1）（2）3有限长序列线性卷积子程序，用于计算两个给定长度（分别为M和N）的序列的卷积，输出序列长度为L=M+N-1。（二）上机实验内容在编制以上各部分程序以后，编制主程序调用各个功能模块实现对信号、系统和系统响应的时域和频域分析，完成以下实验内容。1分析理想采样信号序列的特性产生理想采样信号序列，使：（1）首先选用采样频率为1000Hz，T=1/1000，观察所得理想采样信号的幅频特性，在折叠频率以内和给定的理想幅频特性无明显差异，并作记录。（2）改变采样频率为300Hz，T=1/300，观察所得理想采样信号的幅频特性曲线的变化，并作记录。（3）进一步减小采样频率为200Hz，T=1/200，观察频谱混叠现象是否明显存在，说明原因，并记录此时的幅频特性曲线。2离散信号、系统和系统响应的分析（1）观察信号和系统的时域和频域特性；利用线性卷积求信号通过系统以后的响应。比较系统响应和系统的时域和幅频特性。注意它们之间有无差异，绘出图形。（2）观察信号和系统的时域和幅频特性，利用线性卷积求系统响应。判断响应序列图形及序列非零值长度是否与理论结果一致，说出一种定性判断响应序列图形正确与否的方法（提示：）。利用序列的傅立叶变换数值计算子程序求出，观察响应序列的幅频特性。定性判断结果是否正确。改变信号的矩形宽度，使N=5,重复以上动作，观察变化，记录改变参数前后的差异。（3）将（2）中的信号换为，其中，重复（2）中的实验各步；改变参数再重复（2）中的实验各步；改变参数重复（2）中的实验各步。观察参数的改变对信号与系统响应的时域和幅频特性的影响，绘制相应的图形。3卷积定律的验证。利用式（1-8）将和系统的傅氏变换相乘，直接求得，将得到的幅频特性曲线和实验2（3）中得到的曲线进行比较，观察二者有无差异，验证卷积定律。（三）MATLAB上机内容1阅读本实验指导书中有关MATLAB进行数字信号处理部分，熟悉MATLAB下信号处理的过程和方法。2在MATLAB下重复（二）上机内容的所有要求，将MATLAB的输出结果同自己程序的输出结果进行比较。3改变信号中的衰减因子，先定性估计频谱可能发生的变化，然后观察其频谱的变化，记录结结果，变化是否与你所想的一致？着这说明了什么？4一个LTI系统的冲激响应为，输入序列为，求系统响应和输出信号及其频谱；如果，结果又如何呢？5编写一个程序，将分解为奇偶序列，绘制奇偶序列时域图形并求出它们频谱和，同的频谱进行比较，可以得出什么结论？四思考题1回答上机内容2（2）中的问题。2在分析理想采样信号序列的特性实验中，利用不同采样频率所得到的采样信号序列的傅氏变化频谱，数字频率度量是否相同？它们所对应的模拟频率是否都相同？3在卷积定律的验证实验中，如果选用不同的M值，例如选M=50和M=30，分别作序列的傅氏变换，并求得，所得的结果之间有何差异，为什么？五实验报告要求1在实验报告中简述实验目的和实验原理要点。2总结在上机实验内容中要求比较时域、幅频曲线差异部分内容的结果，定性分析它们正确与否，并简要说明这些结果的含义。3总结实验中的主要结论。4回答思考题。5总结一下你在用MATLAB进行数字信号处理实验项目的时候常用的函数及其功能。6在用MATLAB 处理时和你自己的程序实验时结果是否一致？对不一致的情况进行一个简要的分析。实验二 应用FFT对信号进行频谱分析一实验目的1在理论学习的基础上，通过本次实验，加深对FFT的理解，熟悉FFT算法及其程序的编写。2熟悉应用FFT对典型信号进行频谱分析的方法。3了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。二、实验原理与方法一个连续信号的频谱可以用它的傅立叶变换表示为 (2-1)如果对该信号进行理想采样，可以得到采样序列 (2-2)同样可以对该序列进行Z变换，其中T为采样周期 (2-3)当的时候，我们就得到了序列的傅立叶变换 (2-4)其中为数字角频率，和模拟域频率的关系为 (2-5)式中的是采样频率。上式说明数字角频率是模拟频率对采样速率的归一化。同模拟域的情况相似，数字角频率代表了序列值变化的速率，而序列的傅立叶变换称为序列的频谱。序列的傅立叶变换和对应的采样信号频谱具有下式的对应关系： (2-6)即序列的频谱是采样信号频谱的周期延拓。从式(2-6)可以看出，只要分析采样序列的频谱，就可以得到相应的连续信号的频谱。在各种信号序列中，有限长序列在数字信号处理中占有很重要的地位。无限长的序列也往往可以用有限长序列来逼近。对于有限长的序列我们可以使用离散傅立叶变换（DFT），这一变换可以很好的反映序列的频域特性，并且容易利用快速算法在计算机上实现。当序列的长度为N时，定义DFT为： (2-7)其中，它的反变换定义为： (2-8)令，则有： (2-9)可以得到，是Z平面单位圆上幅角为的点，就是将单位圆进行N等分以后第K个点。所以，X(K)是Z变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。时域采样在满足Nyquist定理时，就不会发生频谱混叠。DFT是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。如同理论课教材所讨论的，在运用DFT进行频谱分析的时候可能有三种误差，即：（1）混叠现象从中可以看出，序列的频谱时采样信号频谱的周期延拓，周期是，因此当采样速率不满足定理Nyquist，经过采样就会发生频谱混叠。这导致采样后的信号序列频谱不能真实的反映原信号的频谱。所以，在利用DFT分析连续信号频谱的时候，必须注意这一问题。避免混叠现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。这告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。（2）泄漏现象 实际中的信号序列往往很长，甚至是无限长。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的DFT来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个矩形窗函数。值得一提的是，泄漏是不能和混叠完全分离开的，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。（3）栅栏效应因为DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度看， 用DFT来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上看到真是的频谱。这样的话就会有一些频谱的峰点或谷点被“栅栏”挡住，不能被我们观察到。减小栅栏效应的一个方法是在源序列的末端补一些零值，从而变动DFT的点数。这种方法的实质是改变了真是频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的一些频谱的峰点或谷点显露出来。注意，这时候每根谱线所对应的频和原来的已经不相同了。从上面的分析过程可以看出，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减小和消除这些误差的影响。FFT并不是DFT不相同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次的分解，使其成为若干小店数DFT的组合，从而减小运算量。常用的FFT是以2为基数的，其长度为。它的运算效率高，程序比较简单，使用也十分的方便。当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。IFFT一般也可以通过FFT程序来完成，比较式和，只要对取共轭，进行FFT运算，然后再去共轭，并乘以因子，就可以完成IFFT。三实验内容及步骤(一) 编制实验用主程序及相应子程序1在实验之前，认真复习DFT和FFT有关的知识，阅读本实验原理和实验附录部分中和本实验有关的子程序，掌握子程序的原理并学习调用方法。2编制信号产生子程序及本实验的频谱分析主程序。实验中需要用到的基本信号包括：（1）高斯序列：（2）衰减正旋序列：（3）三角波序列（4）反三角序列（二）上机实验内容1观察高斯序列的时域和频域特性（1）固定信号中的参数，改变的值，使分别等于2，4，8。观察它们的时域和幅频特性，了解取不同值的时候，对信号时域特性和幅频特性的影响。（2）固定，改变，使分别等于8，13，14，观察参数变化对信号序列时域及幅频特性的影响。注意等于多少时，会发生明显的泄漏现象，混叠现象是否也随之出现？记录实验中观察到的现象，绘制相应的时域序列和幅频特性曲线。2观察衰减正旋序列的时域和幅频特性（1）令并且，检查谱峰出现的位置是否正确，注意频谱的形状，绘制幅频特性曲线。（2）改变，再变化，观察这两种情况下，频谱的形状和谱峰出现的位置，有无混叠和泄漏现象出现？说明出现的原因。3观测三角波序列和反三角波序列的时域和幅频特性（1）用8点FFT分析信号和的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘制两者的序列和幅频特性曲线。（2）在和的末尾补零，用16点FFT分析两个信号的幅频特性，观察幅频特性发生了什么变化？两个信号之间的FFT频谱还有没呀相同之处？这些变化说明了什么？（三）MATLAB 上机内容1在MATLAB下重复（二）中上机实验内容的所以要求，将MATLAB的输出结果同自己程序的输出结果进行比较。2将信号的长度N设为63，用MATLAB中randn（1，N）函数产生一个噪声信号w(n)，计算将这个噪声信号叠加到上以后新信号的频谱，观察发生的变化并记录。3在步骤2的基础上，改变参数和，观察在出现混叠现象和泄漏现象的时候有噪声的信号的频谱有什么变化，是否明显？四思考题1实验中的信号序列和，在单位圆上的Z变化频谱和和会相同吗？如果不同，你能说出哪一个低频分量更多一些吗？为什么？2对一个有限长序列进行离散傅立叶变换（DFT），等价于将该序列周期延拓后进行傅立叶级数展开(DFS)。因为DFS也只是取其中一个周期来运算，所以FFT在一定条件下也可以用以分析周期信号序列。如果实正弦信号，用16点的FFT来作DFS运算，得到的频谱是信号本身的真是谱吗?五实验报告要求1在实验报告中简述实验目的和实验原理要点。2在实验报告中附上在实验过程中记录的各个信号序列的时域和幅频特性曲线，分析所得到的结果图形，说明各个信号的参数变化对其时域和幅频特性的影响。3总结实验中的结论。4回答思考题。实验三 用双线性变换法设计IIR滤波器一实验目的1了解工程上两种最常用的变换方法：脉冲响应不变法和双线性变换法。2掌握双线性变换法设计IIR滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR数字滤波器的计算机程序。3观察用双线性变换法设计的滤波器的频域特性，并与脉冲响应不变法相比较，了解双线性变换法的特点。4熟悉用双线性变换法设计数字Butterworth和Chebyshev 滤波器的全过程。5了解多项式乘积和多项式乘方运算的计算机编程方法。二实验原理与方法从模拟滤波器设计IIR数字滤波器具有四种方法：微分差分变换法、脉冲响应不变法、双线性变换法、z平面变换法。工程上常用的是其中的两种：脉冲响应不变法、双线性变换法。脉冲响应不变法需要经历如下基本步骤：由已知系统传输函数H(S)计算系统冲激响应h(t)；对h(t)等间隔采样得到h(n)=h(n T)；由h(n)获得数字滤波器的系统响应H(Z)。这种方法非常直观，其算法宗旨是保证所设计的IIR滤波器的脉冲响应和模拟滤波器的脉冲响应在采样点上完全一致。而双线性变换法的设计准则是使数字滤波器的频率响应与参考模拟滤波器的频率响应相似。脉冲响应不变法一个重要的特点是频率坐标的变换是线性的（），其确定是有频谱的周期延拓效应，存在频谱混叠的现象。为了克服脉冲响应不变法可能产生的频谱混叠，提出了双线性变换法，它依靠双线性变换式：， ， 其中 ，建立其S平面和Z平面的单值映射关系，数字域频率和模拟域频率的关系是：， （31）由上面的关系式可知，当时，终止在折叠频率处，整个轴单值的对应于单位圆的一周。因此双线性变换法不同于脉冲响应不变法，不存在频谱混叠的问题。从式（31）还可以看出，两者的频率不是线性关系。这种非线性关系使得通带截至频率、过渡带的边缘频率的相对位置都发生了非线性畸变。这种频率的畸变可以通过预畸变来校正。用双线性变换法设计数字滤波器时，一般总是先将数字滤波器的个临界频率经过式（31）的频率预畸变，求得相应参考模拟滤波器的个临界频率，然后设计参考模拟滤波器的传递函数，最后通过双线性变换式求得数字滤波器的传递函数。这样通过双线性变换，正好将这些频率点映射到我们所需要的位置上。参考模拟滤波器的设计，可以按照一般模拟滤波器设计的方法，利用已经成熟的一整套计算公式和大量的归一化设计表格和曲线。这些公式、表格主要是用于归一化低通原型的。通过原型变换，可以完成实际的低通、带通和高通滤波器的设计。在用双线性变换法设计滤波器的过程中，我们也可以通过原型变换，直接求得归一化参考模拟滤波器原型参数，从而使得设计更加简化。理论课教材给出了IIR低通、带通和高通滤波器设计双线性原型变换公式的总结，请参阅。在本实验中，我们只设计两种滤波器（Butterworth和Chebyshev）的设计，相应的这两种参考模拟原型滤波器的设计公式见理论课教材。综上所述，以低通数字滤波器设计为例，可以将双线性变换法设计数字滤波器的步骤归纳如下：1确定数字滤波器的性能指标。这些指标包括：通带、阻带临界频率；通带内的最大衰减；阻带内的最小衰减；采样周期T。2确定相应的数字频率，。3计算经过频率预畸变的相应参考模拟低通原型的频率，4计算低通原型阶数N；计算3dB归一化频率，从而求得低通原型的传递函数。5用表（31）中所列变换公式，代入，求得数字滤波器的传世函数。6分析滤波器频域特性，检查其指标是否满足要求。三实验内容及步骤（一）编制实验用主程序和相应子程序1在实验前复习数字信号处理理论课中有关滤波器设计的知识，认真阅读本实验的原理部分，读懂本指导书附录A4“滤波器有关算法”中列举的有关算法。有理分式代换子程序BSF、多项式乘方展开子程序PNPE和多项式乘积展开子程序YPMP，掌握模拟滤波器传递函数产生子程序BCG、有理分式代换子程序BSF的调用方法。2编制一个用双线性变换法设计IIR数字和滤波器的通用程序。采样周期、通带和阻带临界频率以及相应的衰减等参数在程序运行时输入；根据这些输入参数，计算阶数N、传递函数；输出分子分母系数；绘制幅频特性曲线，绘制点数为50点，()3实验用子程序介绍。附录A4部分中编排的子程序较多，下面简要介绍一下，大家首先要弄清程序的功能、作用，函数的调用参数及意义，再具体了解子程序的编写。（1）模拟滤波器传递函数产生子程序BCG该程序用于产生表中两种逼近方法的传递函数，根据实际需要，产生的传递函数有两种逼近形式的输出：分母以根的形式；分母以多项式的形式。具体见下式：采样Butterworth逼近，还是Chebyshev逼近在子程序运行过程中，由同学们自己选择。调用该子程序的时候，需要输入滤波器的临界频率、和相应衰减、。根据这些条件，子程序可以计算阶数N，产生传递函数。调用格式为：Bcg(AP,AS,WP,WS,&N,*H,*B,&C);其中参数的具体说明如下：AP,AS是浮点数，对应通带衰减、阻带衰减。 WP,WS是浮点数，对应通带临界频率、阻带临界频率。 N是浮点数，存放模拟滤波器的阶数 H是浮点型一维数组，提供分母为多项式形式的传递函数的分母多项式的系数，H(0)=b0，H(1)=b1，H(N)=bN B是浮点型一维复数数组，存放H(S)的分母多项式的根，B(0)=s1，B(1)=s2，B(N-1)=sN C是浮点型数，存放H(S)的分子常数a0，和用于提供分母为根形式的传递函数。（2）有理分式代换子程序BSF该子程序用于双线性变换中的有理分式代换 (3-3)这里只考虑传递函数分子为常数1的情况，这也是通常低通原型的形式。子程序的调用格式为：Bsf(*C,NO,*F1,*F2,NF,*A,*B,&Nob)其中具体参数含义如下：C：浮点型一维数组，存放H(S)分母多项式系数NO:整型数，存放H(S)分母多项式阶次F1:浮点型一维数组，存放变换式的分子多项式系数，对应式(3-3)中F1(z)F2:浮点型一维数组，存放变换式的分母多项式系数，对应式(3-3)中F2(z)NF:整型数，存放变换式s=F1(z)/ F2(z)的阶次A:浮点型一维数组，存放变换结果H(Z)的分子多项式系数B:浮点型一维数组，存放变换结果H(Z)的分母多项式系数Nob:整型数，存放H(Z)的阶次（3）多项式乘方展开子程序PNPE用于多项式的乘方展开，见下式所示：调用格式为：Pnpe(*A,M,N,*B,&MNB)具体参数说明：A:浮点型一维数组，存放多项式的系数数组M: 整型数，存放多项式的阶次N: 整型数，存放乘方次数B: 浮点型一维数组，输出参数，存放展开结果MNB: 整型数，输出展开以后多项式的阶次（4）多项式乘积展开子程序YPMP用于两个多项式的乘积展开，见下式所示：调用格式为：Ypmp(*A,M,*B,N,*C,&MN)具体参数说明：A:浮点型一维数组，存放第一个多项式的系数M:整型数，存放第一个多项式的阶次mB:浮点型一维数组，存放第二个多项式的系数N:整型数，存放第二个多项式的阶次nC:浮点型一维数组，存放结果多项式的系数MN:整型数，存放结果多项式的阶次（二）上机实验内容1采样频率为1Hz，设计一个Chebyshev高通数字滤波器，其中通带临界频率，通带内衰减小于0.8dB，阻带临界频率，阻带内衰减大于20dB。求这个数字滤波器的传递函数H(Z)，输出它的幅频特性曲线，观察其通带衰减和阻带衰减是否满足要求。2采样频率为1Hz，设计一个低通数字滤波器，其中通带临界频率，通带内衰减小于1dB，阻带临界频率，阻带内衰减大于25dB。求这个数字滤波器的传递函数H(Z)，输出它的幅频特性曲线，观察其通带衰减和阻带衰减是否满足要求。3设计一个Butterworth带通数字滤波器，其上、下边带1dB处的通带临界频率分别为20kHz和30kHz，当频率低于15kHz时，衰减要大于40dB，采样周期为10微妙，求这个数字滤波器的传递函数，输出它的幅频特性曲线，观察其通带衰减和阻带衰减是否满足要求。（三）MATLAB 上机内容1在MATLAB下重复（二）上机实验内容的所有要求，将MATLAB的输出结果同自己程序的输出结果进行比较。2在（二）上机实验内容部分只要求用双线性变换法实现IIR滤波器，在MATLAB下将已经实现的滤波器用冲激响应不变法再实现一遍，并将得到的结果进行比较。四思考题1双线性变换法和冲激响应不变法相比较，有哪些优点和缺点？为什么？2双线性变换是一种非线性变换，在实验中你观测到了这种非线性关系吗？应该怎样从哪种数字滤波器的幅频特性曲线中取观察这种非线性关系？3在用MATLAB实现冲激响应不变法设计滤波器的时候，你有没有观察到频谱混叠的现象？请解释。五实验报告要求1在实验报告中简述实验目的和实验原理要点。2记录在上机内容中所设计的滤波器的传递函数H(z)及对应的幅频特性曲线定性分析它们的性能，判断设计是否满足要求。3总结实验中的主要结论。4回答思考题。实验四 用窗函数设计FIR滤波器一实验目的1熟悉FIR滤波器的设计基本方法2掌握用窗函数设计FIR数字滤波器的原理与方法，熟悉相应的计算机高级语言编程。3熟悉线性相位FIR滤波器的幅频特性和相位特性。4了解各种不同窗函数对滤波器性能的影响。二实验原理与方法（一）FIR滤波器的设计在前面的实验中，我们介绍了IIR滤波器的设计方法并实践了其中的双线性变换法，IIR具有许多诱人的特性；但如此同时，也具有一些缺点。例如：若想利用快速傅立叶变换技术进行快速卷积实现滤波器，则要求单位脉冲响应是有限长的。此外，IIR滤波器的优异幅度响应，一般是以相位的非线性为代价的，非线性相位会引起频率色散。FIR滤波器具有严格的相位特性，这对于语音信号处理和数据传输是很重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。本实验中的窗函数法比较简单，可应用现成的窗函数公式，在技术指标要求不高的时候是比较灵活方便的。它是从时域出发，用一个窗函数截取理想的得到，以有限长序列近似理想的；如果从频域出发，用理想的在单位圆上等角度取样得到，根据得到将逼近理想的，这就是频率取样法。（二）窗函数设计法同其它的数字滤波器的设计方法一样，用窗函数设计滤波器也是首先要对滤波器提出性能指标。一般是给定一个理想的频率响应，使所设计的FIR滤波器的频率响应去逼近所要求的理想的滤波器的响应。窗函数法设计的任务在于寻找一个可实现（有限长单位脉冲响应）的传递函数去逼近。我们知道，一个理想的频率响应的傅立叶反变换所得到的理想单位脉冲响应往往是一个无限长序列。对经过适当的加权、截短处理才能得到一个所需要的有限长脉冲响应序列。对应不同的加权、截短，就有不同的窗函数。所要寻找的滤波器脉冲响应就等于理想脉冲响应和窗函数的乘积，即 由此可见，窗函数的形状就决定了滤波器的性质。例如：窗函数的主瓣宽度决定了滤波器的过渡带宽；窗函数的旁瓣大小决定了滤波器的阻带衰减。以下是几种常见的窗函数：1矩形窗2Hannning窗3Hamming窗4Blackman窗5Kaiser窗其中是零阶贝塞尔函数，可以通过改变参数，改变其主瓣宽度和旁瓣大小。在实际设计过程中，上述几种窗函数可以根据对滤波器过渡带宽带和阻带衰减的要求，适当选取窗函数的类型和长度N，以得到比较满意的设计效果。如何根据滤波器长度N的奇偶性，选择的奇偶对称性则是另外一个需要考虑的问题。线性相位实系数FIR滤波器按其N值奇偶和的奇偶对称性，可以分为四种，它们具有不同的幅频和相位特性：1为偶对称，为奇数它的幅度是关于点成偶对称。2为偶对称，为偶数它的幅度是关于点成奇对称，处有零点，所以它不适合用于高通滤波器。3为奇对称，为奇数它的幅度是关于点成奇对称，在处都有零点，所以它不适合用于低通和高通滤波器。4为奇对称，为偶数它的幅度是关于点成奇对称，处有零点，所以它不适合用于低通滤波器。在滤波器设计过程中，只有根据上述四种线性相位滤波器传递函数的性质，合理的选择应采用的种类，构造出的幅频特性和相位特性，才能求得所需要的、具有单位脉冲响应的线性相位FIR滤波器传递函数。窗函数法设计线性相位FIR滤波器可以按照如下步骤进行：（1） （1） 确定数字滤波器的性能要求。确定个临界频率和滤波器单位脉冲响应长度N。（2） （2） 根据性能要求和N值，合理选择单位脉冲响应的奇偶对称性，从而确定理想频率响应的幅频和相位特性。（3） （3） 用傅立叶反变换公式（42），求得理想单位脉冲响应。（4） （4） 选择适当的窗函数w(n)，根据式（43），求得所设计的FIR滤波器单位脉冲响应。（5） （5） 用傅立叶变换求得其频率响应，分析它的幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述过程，直至得到满意的结果。注意：在上述步骤（3）中，从到的反变换要用到式（42），这里的积分运算，在计算机上可取其数值解其中，而，这样，数值解才能较好的逼近解析解。三实验内容及步骤（一）编制实验用主程序及相应子程序1在实验编程之前，认真复习有关FIR滤波器设计的有关知识，尤其是窗函数的有关内容，阅读本次实验指导，熟悉窗函数及四种线性相位FIR滤波器的特性，掌握窗函数设计滤波器的具体步骤。2编制窗函数设计FIR滤波器的主程序及相应子程序。（1）傅立叶反变换数值计算子程序，用于计算设计步骤（3）中的傅立叶反变换，给定，。按照公式求得理想单位脉冲响应，。（2）窗函数产生子程序，用于产生几种常见的窗函数序列。本实验中要求产生的窗函数序列有：矩形窗、Hanning窗、Hamming窗、Blackman窗、Kaiser窗。根据给定的长度N，按照公式（44）到公式（48）生成相应的窗函数序列。在产生Kaiser窗函数的时候需要用到零阶贝塞尔函数，在实验指导书的附录部分提供了零阶贝塞尔函数计算子程序。（3）主程序，在上述子程序的基础上，设计主程序完成线性相位FIR滤波器的窗函数法设计。其中理想滤波器幅频特性的一半（从到区间）频率点上的值，以及滤波器的长度N可以从数据文件或其它形式输入。的另外一半（从到区间）的幅频特性和全部相位特性在程序中根据N的奇偶性和幅频特性的要求，在四种滤波器中选择一种，自动产生。（二）上机实验内容在计算机上调试自己设计好的窗函数法设计FIR线性相位滤波器设计程序。以下是一个例题及其标准答案，供同学们在设计过程中参考，若有不清楚的地方，请大家在翻阅理论课教材。例：用窗函数法设计一个长度N为8的线性相位FIR滤波器，其理想的幅频特性为：分别用矩形窗、Hanning窗、Hamming窗、Blackman窗、Kaiser窗（）设计该滤波器，比较设计结果，下表是参考答案。表（41）程序调试成功以后，请完成以下实验内容：（1） （1） 用Hanning窗设计一个线性相位带通滤波器，其长度N为15，上、下边带截至频率分别为，求，绘制它的幅频和相位曲线，观察它的实际3dB和20dB带宽。如果N为45，重复这一设计，观察幅频和相位特性曲线的变换情况，注意长度N对曲线的影响。（2） （2） 改用矩形窗和Blackman窗，设计步骤（1）中的带通滤波器，观察并记录窗函数对滤波器幅频和相位特性的影响，比较这三种窗函数的特点。（3） （3） 用Kaiser窗设计一个专用的线性相位滤波器。N为40，理想的幅频特性如下图所示：当值分别为4，6，8时，设计相应的滤波器，比较它们的幅频和相位特性，观察并分析值不同的时候对结果有什么影响。（三）MATLAB上机内容1在MATLAB下重复（二）上机实验内容的所有要求，将MATLAB的输出结果同自己程序的输出结果进行比较。2在（二）上机实验内容部分只要求了用窗函数实现FIR滤波器，有兴趣的同在MATLAB下将已经实现的滤波器用频率取样法可再实现一遍，并进行结果比较。四思考题1定性的说明用本实验程序设计的FIR滤波器的3dB截至频率再什么位置？它等于理想频率响应的截至频率吗？2如果没有给定长度N，而是给定了通带边缘截至频率和阻带临界频率，以及相应的衰减，你能根据这些条件用窗函数法设计线性相位FIR低通滤波器吗？3频率取样法和窗函数法各有什么特点？简单说明。五实验报告要求1在实验报告中简述实验目的和实验原理要点。2再试验报告中附上实验过程中记录的h(n)的幅频和相位特性曲线，比较它们的性能，说明滤波器N和窗函数对滤波器性能的影响。3总结实验中窗函数法设计FIR的特点，归纳设计中的主要公式。4回答思考题。附录A C语言实现数字信号处理算法 附录A1 BC下复数类型的实现 1、利用BC提供的复数支持 /BC中使用复数类型使用示例(ComplexUse.Cpp文件) #include #include int main(void) double x = 3.1, y = 4.2; complex z = complex(x,y); cout z = z n; cout and imaginary real part = imag(z) n; cout z has complex conjugate = conj(z) n; return 0; 2、定义复数类，填写相应成员函数 /C中的复数类型调用时可能不是非常好用，可自己定义复数类(ComplexUse.Cpp文件) class Complex public: Complex() Complex( float re, float im ); float r()return real; float i()return imag; float mod()return sqrt(real*real+imag*imag); Complex operator+( Complex &other ); Complex operator-( Complex &other ); Complex operator*( Complex &other ); Complex operator/( Complex &other ); private: float real, imag; ;/ Operator overloaded using a member function Complex:Complex(float re,float im) real=re; imag=im; ; Complex Complex:operator+( Complex &other ) return Complex( real + other.real, imag + other.imag ); ; Complex Complex:operator-( Complex &other ) return Complex( real - other.real, imag - other.imag ); ; Complex Complex:operator*( Complex &other ) float x,y; x=real*other.real-imag*other.imag; y=real*other.imag+imag*other.real; return Complex( x,y ); ; Complex Complex:operator/( Complex &other ) float x,y,l; l=other.real*other.real+other.imag*other.imag; x=real*other.real+imag*other.imag; y=other.real*imag-real*other.imag; x=x/l; y=y/l; return Complex(x,y); ; 附录A2 BC下的绘图 1、通用绘图函数的有关说明 （1）函数输入参数说明 left : 绘图区的左上角横坐标 top : 绘图区的左上角纵坐标 right : 绘图区的右下角横坐标 bottom : 绘图区的右下角纵坐标 f : 需要绘制图形的数组 length : f数组的长度 （2）函数调用示例 /在以(10,5)为左上角坐标，(200,240)为右下角坐标的区域内绘制数组x的图形，数组长度为10 Plot(10,5,200,240,x,10); （3）BC图形程序需要“BC安装目录bgiegavga.bgi”文件支持 2、通用绘图函数 void Plot(int left,int top,int right,int bottom,double *f,int length) int n; int Left,Top,Right,Bottom,Mid; double x_unit_length,y_unit_length,max=-1.0; for(n=begin;n=end;n+) /begin和end表示数组下标的起始值和结束值 if(maxfabs(fn) max=fabs(fn); Left=left+5;Top=top+5; Right=right-5;Bottom=bottom-5; Mid=(Bottom+Top)/2; x_unit_length=(Right-Left)/(N+2); y_unit_length=(Bottom-Mid)/max; int gdriver = DETECT, gmode; initgraph(&gdriver, &gmode, ); rectangle(left,top,right,bottom); setcolor(YELLOW); setcolor(GREEN); int x; x=Left+x_unit_length/2; n=0; while(nlength) line(x,Mid,x,Mid-fn*y_unit_length); x+=x_unit_length; n+; setcolor(LIGHTGRAY); getch(); closegraph(); 附录A3 傅立叶变换有关算法 1、傅立叶变换数值计算函数 /* 输入参数： f : 需要进行傅立叶变换的数值序列 N : 输入序列的长度 M : 频谱分析的频点数 result : 傅立叶变换的结果（复数序列） 调用示例： M=100; DSFT(f,0,N,M,F); */ void DSFT(double *f,int N,int M,complex *result) int k,n; double omega,delta_omega; delta_omega=2*M_PI/M; omega=0.0; for(k=0;kM;k+) resultk=complex(0.0,0.0); for(n=0;nN;n+) resultk+=fn*exp(complex(0,-1)*omega*n); omega+=delta_omega; 2、卷积计算函数 /* 参数说明 f1 : 序列一； f1_begin : 序列一开始下标； f1_end : 序列一结束下标； f2 : 序列二； f2_begin : 序列二开始下标； f2_end : 序列二结束下标 f : 存放结果序列； begin: 存放结果序列开始下标； end : 存放结果序列结束下标 调用示例： JuanJi(f1,0,10,f2,0,3,f3,&begin,&end); */ void JuanJi(double *f1,int f1_begin,int f1_end,double *f2,int f2_begin,int f2_end,double *f,int *f_begin,int *f_end) int m,n,i,begin,end; *f_begin=(f1_beginf2_begin)?f1_begin:f2_begin; *f_end=*f_begin+f1_end-f1_begin+f2_end-f2_begin; for(i=0,n=*f_begin;nf1_begin)?(n-f2_end):f1_begin; end=(n-f2_beginf1_end)?(n-f2_begin):f1_end; fi=0; for(m=begin;m=end;m+) fi+=f1m-f1_begin*f2n-m-f2_begin; 3、FFT变换函数 /* 参数说明： xxx : 需要进行变换的数组； N : 长度； FFtorIFFT : 进行正变换或反变换的指示参数 调用示例： for(i=1;i=256;i+=2) xxxi=1.0; /xxx1,3,5,7. is real part xxxi+1=0.0; /xxx2,4,6,8. is image part FFT(xxx,128,-1); */ #define SWAP(x,y) double temp=(x);(x)=(y);(y)=temp; void FFT(double xxx,int N,int FFTorIFFT)/*FFT 1 IFFT -1*/ double temp,tempr,tempi; double wtemp,wr,wpr,wpi,wi,theta; int n,mmax,m,q,istep,p; n=N1; q=1; for(p=1;pp) SWAP(xxxq,xxxp); SWAP(xxxq+1,xxxp+1); m=n1; while(m=2&qm)q-=m;m=1; q+=m; mmax=2; while(nmmax) istep=2*mmax; theta=2.0*M_PI/(FFTorIFFT*mmax); wtemp=sin(double)(0.5*theta); wpr=-2.0*wtemp*wtemp; wpi=sin(double)theta); for(wr=1.0,wi=0.0,m=1;mmmax;m+=2) for(p=m;p=n;p+=istep) q=p+mmax; tempr=wr*xxxq-wi*xxxq+1; tempi=wr*xxxq+1+wi*xxxq; xxxq=xxxp-tempr; xxxq+1=xxxp+1-tempi; xxxp+=tempr; xxxp+1+=tempi; wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; mmax=istep; 附录A4 滤波器有关算法 1、多项式的乘方展开函数 /* 参数说明： A : 一维数组，为多项式的系数数组 M : 多项式的阶次 N : 乘方次数 B : 一维数组，输出参数，存放展开结果 MNB : 输出展开以后多项式的阶次 */ int Pnpe(double *A,int M,int N,double *B,int &MNB) double C20; int i,NK,k,j; MNB=M*N; for (i=1;i20;i+) Ci=0.0; Bi=0.0; if(N=0)B1=1.0;return 0; for(i=1;i=M+1;i+) Bi=Ai; if(N=1) return 0; NK=M+1; for(i=2;i=N;i+) for (j=1;j=M+1;j+) for(k=1;k=NK;k+) Ck+j-1=Ck+j-1+Aj*Bk; NK=M+NK; for(k=1;k=NK;k+)Bk=Ck; Ck=0.0; return 0; 2、多项式乘积展开函数 /* 参数说明： A : 一维数组，存放第一个多项式的系数 M : 第一个多项式的阶次 B : 一维数组，存放第二个多项式的系数 N : 第二个多项式的阶次 C : 一维数组，存放结果多项式的系数 MN : 结果多项式的阶次 */ int Ypmp(double *A,int M,double *B,int N,double *C,int &MN) int i,j; MN=M+N; for (i=1;i=20;i+) Ci=0.0; for(i=1;i=M+1;i+) for (j=1;j=N+1;j+) Ci+j-1=Ci+j-1+Ai*Bj; return 0; 3、有理分式替换函数 /* 参数说明： C : 一维数组，存放H(s)分母多项式系数 NO : H(s)分母多项式阶次 F1 : 一维数组，存放变换式的分子多项式系数 F2 : 一维数组，存放变换式的分母多项式系数 NF : 变换式的阶次 A : 一维数组，变换结果的分子多项式系数 B : 一维数组，变换结果的分母多项式系数 Nob : H(z)的阶次 */ int Bsf(double *C,int NO,double *F1,double *F2,int NF,double *A,double *B,int &Nob) double D20,E20,G20; int I,j,ii,NND,NNE,NG; Pnpe(F2,NF,NO,A,Nob); for (I=1;I=20;I+) BI=0.0; for (I=1;I=NO+1;I+) ii=I-1; Pnpe(F1,NF,ii,D,NND); Pnpe(F2,NF,NO-ii,E,NNE); Ypmp(D,NND,E,NNE,G,NG); for(j=1;j=NG+1;j+) Bj=Bj+CI*Gj; return 0; 4、Butterworth滤波器阶数计算函数 int BN (double alpha_p, double alpha_s, double omiga_p, double omiga_s) double t1,t2,t3,t4; t1 = pow(10.0,0.1*alpha_p)-1.0; t2 = pow(10.0,0.1*alpha_s)-1.0; t3 = log10(omiga_p/omiga_s); t4 = log10(sqrt(t1/t2)/t3; return (int)ce