数字信号处理实验讲义.doc

实验一 信号、系统与系统响应南方医科大学数字信号处理实验讲义林霖 杨丰 编著生物医学工程学院 电子技术系二00七年 九月45目录实验一 信号、系统与系统响应1实验二 离散信号的DTFT和DFT7实验三 FFT算法的应用13实验四 脉冲响应不变法和双线性变换法设计IIR数字滤波器19实验五 用窗函数法设计FIR数字滤波器26实验六 信号滤波34附录二 语音信号的基础知识37附录一 MATLAB简介40参考文献42实验一 信号、系统与系统响应一、 实验目的：1. 熟悉数字信号处理中的基本信号，了解信号的时域特性和频域特性，学会利用Matlab仿真软件产生这些基本信号。2. 掌握离散信号的线性卷积，利用卷积方法观察分析系统的时域特性。3. 掌握离散LTI系统特性，观察离散信号通过LTI系统的输出信号，分析输入信号和输出信号之间的差异，加深LTI系统对离散信号的影响的认识。4. 掌握数字滤波器的基本特性，学会根据系统要求设计简单的数字滤波器的方法。二、 实验原理：1. 信号处理中的基本信号Matlab提供了许多函数用于产生常用的基本信号：如阶跃信号、脉冲信号、指数信号、正弦信号和周期方波信号等等。这些基本信号是信号处理的基础。1) 单位脉冲序列单位脉冲序列的定义是： （1.1）位移了n个样本点的单位脉冲序列为： （1.2）2) 单位阶跃序列单位阶跃信号定义是： （1.3）3) 矩形序列长度为N的矩形序列信号是： （1.4）4) 指数序列指数序列定义为： （1.5）右边指数序列是： （1.6）Matlab在信号产生方面有着非常强大的功能，许多函数都可以用来产生这些基本信号，下面列举一些实验中能用到的Matlab函数，函数的具体用法参考Matlab的帮助文件。 plot(x,y)：绘制二维曲线图形，x轴和y轴为线形坐标。用来表示连续的信号曲线。 stem(x,y)：绘制二维杆状图形，x轴和y轴为线形坐标。用来表示离散的信号点。 zeros(m,n)：产生元素全为0的矩阵，其中m为行数，n为列数。 ones(m,n)：产生元素全为1的矩阵，其中m为行数，n为列数。 rectpulse(t,w)：产生高度为1，宽度为w，关于t=0对称的矩形脉冲信号。 square(w0,duty)：产生基本频率为w0，占空比为duty的周期方波信号。 randn(m,n)：产生一个均值为0，方差为1的随机矩阵，其中m为行数，n为列数。2. 信号的基本运算信号的运算是数字信号处理的重要内容。对输入信号按指定的算法进行运算以提取有用的信息，而这些算法可以分解为若干基本运算。离散信号的基本运算主要包括：信号的加、乘、移位、翻转、卷积、相关等，其中，离散信号的卷积运算是我们学习的重点，也是难点。离散信号的线性卷积(convolution)定义是： （1.7）Matlab提供了计算两个离散序列卷积的函数conv，它的调用方式是y=conv(x,h)其中调用参数中x,h为卷积运算所需的两个序列，返回值y是卷积结果。值得注意的是：y序列的起始点是x序列起始点和h序列起始点的和，而y序列的终止点是x序列终止点和h序列终止点的和。另外，conv函数也可以用来计算多项式乘法，具体的用法参考Matlab的帮助文件。3. 离散LTI系统离散LTI系统的输入和输出关系是 （1.8）其中xn，yn，hn分别是系统的输入、输出和单位脉冲响应。系统的频率响应是： （1.9） （1.10）现在重点研究一下系统的响应函数，一般来说，是一个复函数，我们可以把写成模和相位的表达式： （1.11）其中表示的是频率响应函数的模，一般被称作系统的增益。输出信号中不同的频率分量的模等于输入信号中不同的频率分量的模乘以对应频率的系统响应函数的模。即有： （1.12）而则表示的是频率响应函数的相位，一般被称作系统的相移。系统的相移可以改变输入信号中各分量之间的相对相位关系，这样即使系统的增益对所有频率都为常数的情况下，也有可能在输入的时频特性上产生很大的变化，也就是相位失真。（详细的论述参考Oppenheim的Sighals and Systems中的相关内容）。这里给出一个不会引起相位失真的条件：系统的相位在通带范围内尽可能的呈现线性特征。具有线性相位特性的系统有一个特别简单的意义，这就是时移，事实上，相位特性函数的斜率就是时移的大小，而这个斜率的大小就是所谓的群时延。Matlab中提供了计算系统响应的相关函数：用下面的线性常系数差分方程描述的离散LTI系统 （1.13）y=filter(b,a,x)：该函数可以方便地计算出差分方程的零状态响应。其中x表示的是输入信号，y表示通过系统后的输出信号，输出信号的长度和输入信号相同，而，。如果系统函数用的有理多项式来表征，即有 （1.14）则可以用freqz函数来计算系统函数的DTFT抽样值，它的调用形式是：H=freqz(b,a,w)，其中b和a分别代表上式中分子多项式和分母多项式系数的向量，即，。函数中w为抽样的频率点，调用函数时，w的长度至少为2，返回值H就是DTFT在抽样点w上的值，一般是复数。 abs(X)： 求X的幅度 angle(X)： 求X的相位 real(X)： 求X的实部 imag(X)： 求X的虚部4. 简单数字滤波器的设计设计简单的数字滤波器有很多种方法，这里主要讲述一下FIR滤波器。FIR数字滤波器的系统函数可以表征为： （1.15）其中是的次多项式， 如果我们需要设计一个低通滤波器，输入信号的低频分量的角频率是，高频分量的角频率是，为了简单考虑，那么可以假设FIR数字滤波器的长度是2，则有 （1.16）那么系统的差分方程就是： （1.17）系统的频率响应是： （1.18）所以滤波器的群延迟为：为了滤除高频分量，保留低频分量，有 （1.19）解方程得到： （1.20）三、实验内容1. 复习教材第一章的主要内容，认真阅读本讲义中对实验原理的阐述，熟悉matlab下数字信号处理的过程和方法。2. 编写matlab程序，产生下列信号：1) xk=3,2,-1,0,2,-3,0,-1，k=-2:5。2) ，其中=1，采样频率。3) ，其中是均值为0，方差为1的随机信号。3. 已知两个有限长序列：xk=1,2,1,1,0,-3：k=0:5，hk=1,-1,1：k=-1:1。计算离散卷积和，并在图形中表示出信号yk。4. 已知输入信号，其中，是均值为0，方差为1的随机信号；系统函数的差分表达式是，其中。1) 用Matlab画出系统函数的幅频响应和相频响应。2) 画出输入信号和输出信号的曲线图，比较有什么不同。3) 如果输入信号，画出输入信号和输出信号的曲线图，和上一小题比较，看看有什么不同。5. 一个长度为5的FIR滤波器具有对称的脉冲响应且。若输入为三个角频率分别为，的余弦序列的和，写出只能使输入信号中频部分通过的单位脉冲冲激响应，并用Matlab验证系统的滤波效果。四、思考题1. 回答实验内容第3题的问题。2. 实验内容第3题中，输出信号相比于输入信号变得平滑了，请分析原因。3. 实验内容第3题中，输入信号和输出信号会有一个相位的移动，这是为什么？分析原因。4. 总结在实验中使用的主要函数及其功能。实验二 离散信号的DTFT和DFT实验二 离散信号的DTFT和DFT一、实验目的1. 复习离散信号的DTFT和DFT的基本原理，加深对DTFT和DFT的基本性质的理解，学会离散信号频谱分析方法。2. 掌握循环卷积的基本概念，验证循环卷积定理，并会利用DFT计算线性卷积。3. 加深理解连续时间信号离散化过程中的数学概念和物理概念，掌握时域抽样定理，频域抽样定理的基本内容，掌握由抽样序列重建原连续信号的基本原理与实验方法。4. 利用DFT分析连续非周期信号的频谱，理解在引入窗函数的近似分析过程中出现的混叠现象、泄露现象和栅栏现象。二、实验原理1. 有限长序列的离散傅里叶变换(DTF)对于离散非周期信号，其离散时间傅里叶变换是周期为的数字角频率的连续函数。根据频域抽样定理，如果信号为有限长的序列，则可以表示为有限项的虚指数信号线性组合，即有限长的序列的傅里叶表示(DFT)为: （2.1）其中, （2.2）从DFT与DFS之间关系可知，和分别是周期序列和在主值区间上的有限序列，也就是说，当计算周期序列在区间上的数值时，只需要序列在区间上的数值；同样，当计算周期序列在区间上的数值时，只需要序列在区间上的数值。通过对信号在时域进行离散化或周期化等处理，可以建立有限长序列与四种信号在频域的关系，从而实现利用有限长序列的离散傅里叶变换分析其他频谱。2. 循环卷积的基本原理和应用设和均为长度为的有限长序列，两序列点循环卷积定义为： （2.3）两个长度为的有限序列的循环卷积就等于该两个序列周期化后的周期卷积（主值区间上的值）。设和均为长度为的有限长序列，且有 （2.4）则DFT的时域卷积定理为： （2.5）类似的，DFT的频域卷积定理为： （2.6）由于点的和都存在成熟的快速算法，因此利用来实现序列的线性卷积，可以极大地提高线性卷积的效率。长度为的序列与长度为的序列的线性卷积，可以通过两序列的来计算，流程图如下图1所示： 图1 DFT计算线性卷积原理图如果线性卷积的序列中有一个是无穷长的序列或者比另一个长很多的有限序列，可以利用分段卷积的方法处理，主要有两种方法：重叠相加法和重叠保留法。这里讨论重叠相加法。设序列的长度为，为无限长序列。将均匀分段，每段长度取， （2.7）于是有 （2.8）上式说明了，计算与的线性卷积时，可先进行分段线性卷积，然后把分段卷积结果叠加起来即可。Matlab中没有直接调用的循环卷积函数，下面的代码以函数的形式实现了两个序列的循环卷积运算。function fn=circonvt(x1,x2,N)% circonvt函数实现输入序列x1和x2的循环卷积,fn为输出序列%N 为循环卷积长度%实现方法:fn=sum(x1(m)*x2(n-m) mod N)if (length(x1)N | length(x2)N) %判断输入信号的长度 error(N的长度必须大于输入数据的长度);endx1=x1,zeros(1,N-length(x1); %给序列补零x2=x2,zeros(1,N-length(x2); %给序列补零m=0:N-1;x=zeros(N,N);for n=0:N-1 x(:,n+1)=x2(mod(n-m),N)+1); %构造循环矩阵end;fn=x1*x; %循环计算卷积3. 时域抽样与频域抽样计算机等离散系统分析处理连续时间信号时必须对信号进行离散化处理，时域抽样解决了把连续信号转变成适于数字系统处理的离散信号的问题。对连续时间信号以抽样间隔抽样，得到的离散时间信号可以表示为： （2.9）设的频谱为，的频谱为，为抽样角频率，为频率。因为数字角频率与模拟角频率的关系为，有： （2.10）由上式看出，要使频域周期化的过程中不发生频谱混叠，则，这就是时域抽样定理。频域抽样与时域抽样形成对偶关系，是将离散化以便于计算机进行数值计算。在内对进行点均匀抽样，导致时域序列以点为周期进行周期延拓，即： （2.11）当不够大时，将导致时域发生混叠。频域抽样定理给出了频域抽样过程中时域不发生混叠的约束条件：若序列的长度为，则应该有。信号的重建是信号抽样的逆过程，即将离散信号转换为连续时间信号。这需要在样点之间进行内插，记内插函数为，得到的重建信号为，则有： （2.12）4. 利用分析连续非周期信号的频谱通过建立的离散傅里叶变换与连续非周期信号的傅里叶变换之间的关系，可以利用对连续非周期信号频谱进行近似分析，但是，这个近似分析过程中一般会引入三种现象，即混叠现象、泄漏现象和栅栏现象。1) 混叠现象对连续信号进行时域采样，相当于把周期化为，若连续信号为带限信号，采样只要满足时域抽样定理，在中可以得到完整的。这样可以通过频域抽样由来计算。如果连续信号不是带限信号，或者不满足时域抽样定理，在连续信号离散化时，就会出现信号频谱混叠。解决连续信号离散化过程中的频谱混叠主要有两种方法：对于带限连续信号，只要提高抽样频率使之满足时域抽样定理；对于非带限连续信号，可根据实际情况对其进行低通滤波，使之成为带限信号。2) 泄漏现象如果连续信号在时域无限长，则离散化后的序列也无限长，无法适用分析，需要对其进行加窗截短使之成为有限长序列，即： （2.13）加窗处理对信号频谱分析主要有两个方面影响：一是计算出的频谱中有多余的高频分量，这是由窗函数突然截断信号而引起的，这个现象称为频率泄漏；二是谱线变成了具有一定宽度的谱峰，从而降低了频率的分辨率。3) 栅栏现象通过对连续信号抽样加窗得到有限长序列，利用点计算其频谱，由于是离散序列，因而无法反映抽样点之间的全部细节，称之为栅栏现象。栅栏现象是利用分析连续信号频谱过程中无法克服的现象，为了改善栅栏现象，在序列补零是常用的方法。可是，尽管对序列补零似乎可以提高信号的频率分辨率，但这里的频率分辨率只是频谱显示的分辨率。换句话说，如果连续信号在离散化和时域加窗过程中，由于混叠或者泄漏等过程已经造成信号频谱中信息的失真，则无论怎么补零也无法再恢复以损失的信息。5. fft函数Matlab中提供了fft函数来实现DFT的快速算法，其调用格式是： X=fft(x)：基2时间抽取FFT算法，用于计算序列x的傅里叶变换。当x的长度为2的整数次幂或者x全为实数时，计算的时间会大大缩短。 X=fft(x,n)：补零或者截短的n点傅里叶变化。 x=ifft(X)和x=ifft(X,n)是相应的傅里叶反变换。 fftshift(x)：将fft计算输出的零频移到输出的中心。三、实验内容1. 已知序列， (1) 编写Matlab程序，利用freqz函数，计算给序列的表达式，并画出是的曲线。(2) 编写Matlab程序，利用fft函数，计算时，序列的在的抽样值。利用hold函数，将抽样点画在地曲线上。2. 已知序列xk=1,-2,3,1；k=0,1,2,3，hk=3,0,1,-1,2,1；k=0,1,2,3,4,5(1) 编写循环卷积函数circonvt(x1,x2,N)(2) 计算循环卷积。(3) 验证循环卷积定理：。3. 已知序列xk=1,1,1；k=0,1,2，对其频谱进行抽样，分别取N=2，3，10，观察频域抽样造成的混叠现象。4. 已知已离散序列为，。(1) 如果抽样频率为，则，取，其中，。对做64点的，画出此时信号的谱。(2) 如果（1）中显示的谱不能分辨两个谱峰，是否可对（1）中的64点信号补零而分辨出两个谱峰？请通过编程进行验证。(3) 如果增加采样点数，取，那么对做255点的，画出此时信号的频谱。四、思考题1. 在实验内容第一题中，对比分析结果，并解释原因。2. 总结线形卷积定理和循环卷积定理的适用范围，指出那些信号的傅里叶变换适用线形卷积定理，那些信号的傅里叶变换适用循环卷积定理。 3. 解释实验内容第三题的混叠现象。4. 在实验内容第四题中，解释（2）题的实验结果。5. 在实验内容第四题中，（3）题的频谱图是否能分辨出两个谱峰，并对比（1）的结果，解释这种现象的原因。实验三 FFT算法的应用实验三 FFT算法的应用一、实验目的1. 进一步加深DFT算法原理和基本性质的理解，FFT是DFT的一种快速算法，FFT的运算必然满足DFT的基本性质。2. 掌握FFT算法原理，熟悉FFT的算法流程。3. 学习FFT算法的应用。二、实验原理1. 离散信号的傅里叶变换我们先来复习一下离散信号的傅里叶变换，其变换公式有： （3.1）有限长序列的离散傅里叶变换也可以用矩阵形式表示，于是有： (3.2）其中，向量由频域序列的个系数构成，即 （3.3）向量由时域序列的个样本值构成，即 是的矩阵，表示为 （3.4）由定义同样可以得到的矩阵表示为 (3.5）如果利用上面的式子直接计算点序列的，对于每一个频率分量，需要计算次复数乘法，次复数加法。因此，计算个不同频率分量，共需要次复数乘法，次复数加法。的直接计算与的直接计算具有相同的运算量。显然，随着的增大，和的运算量将会急剧的增加。由上面的表述可以看出，的非常适合采用矩阵运算的方式来实现，这种特点为我们利用matlab编写算法程序提供了很大的方便。2. 离散傅里叶变换快速算法在信号处理中得到广泛应用，一个非常重要的原因就是其存在高效算法，即快速傅里叶变换。是计算的各种快速算法的统称。1) 基2时间抽取算法对长度为的序列进行时间抽取，即将其分解为两个长度为点的序列。 （3.6）由的定义可得： （3.7）令 （3.8） （3.9）而由于旋转因子存在对称性，则综合得出： （3.10）由上式可知，将序列按奇偶分解得到两个子序列和，就可以由两个子序列和对应的合成序列的。这种将时域长序列逐次按奇偶分解为两个短序列的称为基2时间抽取。基2时间抽取算法总共需要次复数乘法，而需要的复数加法次数为。根据基2时间抽取算法原理，可以总结出相应的编程思想和程序框图，如下图3-1：图3-1 FFT算法框图2) 基2频率抽取算法基2频率抽取的算法和基2时间抽取的算法原理类似，也是将长序列逐次分解为两个短序列，最后由短序列的合成长序列的。这里不推导算法的过程，直接给出最后表达式设长序列的长度为，分解的两个短序列为 （3.11）它的为 （3.12）3. 算法的应用1) 利用点复序列的同时计算两个点实序列的设和都是点实序列，和分别表示它们对应的点。由点序列和构造一个点复序列，有 （3.13）对复序列进行一次点得： （3.14）根据的线性特征得 （3.15）这样就通过一个点复序列的运算同时得到两个点实序列的2) 利用点复序列的算法计算点实序列的设是一个长度为的实序列，是其点的。定义两个长度为的实序列和为： （3.16）由（3.10）得到 （3.17）3) 利用快速计算利用与定义的对称性，直接由计算 （3.18）上式表明，由通过计算，可以转变为由通过来计算，只需要将的运算结果再取一次共轭，并乘以，就可得到的值。三、实验内容1. 复习直接计算的有关内容，根据矩阵，编写matlab函数程序，以实现的直接运算。2. 已知离散序列1) 利用上题中编写的函数，求序列的256点的离散傅里叶变换。并统计变换所用的时间2) 利用matlab提供的函数，求序列的256点的快速离散傅里叶变换，并统计变换所用的时间。附：统计时间可以利用函数：cputime。调用例子如下： t=cputime; your_operation; runtime=cputime-t；3. 假设通用计算机的速度为平均每次复数乘法需要5um，每次复数加法需要1um，用来计算1024点，问直接计算需要多少时间。如果用计算呢？4. 按照下面的算法： 编写程序，其中部分不用写出清单，可调用子程序。四、思考题1. 已知复序列的8点为V0=1-3j，V1=-2+4j，V2=3+7j，V3=-4-5j，V4=2+5j，V5=-1-2j，V6=4-8j，V7=6j。不计算Vm的，试确定序列和的8点和。2. 在实验内容第3题中，照题目的计算，用进行快速卷积对信号进行处理时，估算可实现实时处理的信号最高频率。3. （选做）复习算法原理和流程图，根据其编程思想和程序框图，试编写的源程序，实现离散序列的快速傅里叶变换。实验四 脉冲响应不变法和双线性变换法设计IIR数字滤波器实验四 脉冲响应不变法和双线性变换法设计IIR数字滤波器一、实验目的1. 加深对数字滤波器的常用指标和设计过程的理解。2. 熟悉用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理与方法。3. 了解多项式乘积和多项式乘方的计算机编程方法。4. 掌握数字滤波器的计算机仿真方法。二、实验原理IIR数字滤波器一般为线性时不变的因果离散系统，阶IIR数字滤波器的系数函数可以表示为的有理多项式，即 （4.1）系数中至少有一个非零。对于因果IIR数字滤波器，满足。1. 滤波器设计中涉及到的技术指标和：分别是数字和模拟滤波器通带截止频率（passband cutoff frequency）和：分别是数字和模拟滤波器阻带截止频率（stopband cutoff frequency）：通带波纹峰值（peak passband ripple）：阻带波纹峰值（peak stopband ripple）：通带的最大衰减（passband attenuation），定义为：阻带的最大衰减（stopband attenuation），定义为：滤波器的增益（gain）响应，定义为2. 数字滤波器设计的基本步骤：a) 将数字滤波器技术指标转换为对应的模拟滤波器技术指标b) 设计满足技术指标的模拟滤波器，c) 将设计出来的模拟滤波器转换为满足技术指标的数字滤波器3. 模拟低通滤波器的设计常用的IIR数字滤波器设计是以模拟滤波器设计为基础的，本实验中，我们只涉及到Butterworth（巴特沃思）和Chebyshev（切比雪夫）两种滤波器的设计，相应的两种参考模拟原形滤波器的设计公式如下：a) Butterworth（巴特沃思）低通滤波器首先由滤波器的设计指标确定滤波器的阶数N （4.2）然后确定滤波器参数 （4.3）最后确定滤波器的系统函数如果N为偶数，则有 （4.4）其中，如果N为奇数，则有 （4.5）其中，b) Chebyshev（切比雪夫）低通滤波器首先由滤波器的设计指标确定滤波器的阶数N （4.6）然后确定滤波器的系统函数先确定几个参数的值 （4.7） （4.8） （4.9）如果N为偶数，则有 （4.10）如果N为奇数，则有 （4.11）c) Chebyshev（切比雪夫）低通滤波器首先由滤波器的设计指标确定滤波器的阶数N （4.12）然后确定滤波器的系统函数先确定几个参数的值 （4.13） （4.14）如果N为偶数，则有 （4.15）如果N为奇数，则有 （4.16）4. 模拟域频率变换在实际应用中，经常需要高通、带通、带阻等其他类型的滤波器。这些滤波器的设计可以通过频率变换的方法转换为低通滤波器来设计，在低通滤波器设计完成以后，再通过频率变换将其转换为所需类型的滤波器。表4.1为变换公式与参数表表4.1 变换公式及参数表变换公式参数注备低通到低通是一正的参数低通到高通是一正的参数低通到带通低通到带阻5. 脉冲响应不变法实现数字滤波器的设计脉冲响应不变法是将模拟滤波器转换为数字滤波器的方法。它通过对模拟滤波器的单位脉冲等间隔抽样来获取数字滤波器的单位脉冲响应。利用脉冲响应不变法实现数字滤波器的设计，滤波器的频率响应都存在一定程度的混叠。对模拟低通和带通滤波器，可通过减小抽样间隔或者增大模拟滤波器在阻带的衰减，使混叠引起的误差在可接受的范围内。对模拟高通和带阻滤波器，由于存在严重的混叠，所以不能用脉冲响应不变法将模拟高通和带阻滤波器转换为数字滤波器。用脉冲响应不变法设计数字滤波器的设计步骤为：1) 利用模拟频率和数字频率的关系将数字滤波器的频率指标转换为模拟滤波器的频率指标；2) 根据模拟滤波器的技术指标设计模拟滤波器的；3) 用脉冲响应不变法，将模拟滤波器转换为数字滤波器。将乘上常数因子，然后对做部分分式展开，对展开中的每一项作变换，从而获得数字滤波器，变换公式为： （4.17）6. 双线性变换法实现数字滤波器的设计双线性变换是一种应用非常广泛的将模拟系统变换成为数字系统的方法。双线性变换的基本思想是利用数值积分将模拟系统变换成为数字系统。双线性变换法只适合设计幅度响应为分段常数的数字滤波器。用双线性变换法设计数字滤波器的设计步骤为：1) 由模拟频率和数字频率的关系将数字滤波器的频率指标转换为模拟滤波器的频率指标；2) 根据模拟滤波器的技术指标设计模拟滤波器的；3) 用双线性变换法，将模拟滤波器转换为数字滤波器。变换公式有： （4.18）7. 常用的IIR滤波器设计的Matlab函数Matlab信号处理工具箱中提供了IIR滤波器设计的函数。常用的函数有：1) Butterworth模拟低通滤波器 N,wc=buttord(wp,ws,Ap,As,s) num,den=butter(N,wc,s)2) Chebyshev型模拟低通滤波器 N,wc=cheb1ord(wp,ws,Ap,As,s) num,den=cheby1(N,Ap,wc,s)3) Chebyshev型模拟低通滤波器 N,wc=cheb2ord(wp,ws,Ap,As,s) num,den=cheby2(N,Ap,wc,s)4) 椭圆模拟低通滤波器 N,wc=ellipord(wp,ws,Ap,As,s) num,den=ellip(N,Ap,As,wc,s)上面的几个滤波器设计函数的参数基本含义是一致的：输入参数wp和ws(rad/s)表示滤波器的通带和阻带截止频率，Ap和As(dB)表示滤波器的通带和阻带衰减，s表示所设计的是模拟滤波器。函数返回参数N表示滤波器的阶数，wc(rad/s)表示滤波器的3dB截止频率。而数组num和dem则表示滤波器系统函数的分子多项式和分母多项式。5) Butterworth数字低通滤波器 N,Wc=buttord(Wp,Ws,Ap,As) num,den=butter(N,Wc)6) Chebyshev型数字低通滤波器 N,Wc=cheb1ord(Wp,Ws,Ap,As) num,den=cheby1(N,Ap,Wc)7) Chebyshev型数字低通滤波器 N,Wc=cheb2ord(Wp,Ws,Ap,As) num,den=cheby2(N,Ap,Wc)8) 椭圆数字低通滤波器 N,Wc=ellipord(Wp,Ws,Ap,As) num,den=ellip(N,Ap,As,Wc)上面的几个滤波器设计函数的参数基本含义是一致的：输入参数Wp和Ws表示数字低通滤波器的归一化通带和阻带截止频率，Ap和As(dB)表示滤波器的通带和阻带衰减。函数返回参数N表示滤波器的阶数，Wc表示滤波器的3dB截止频率。而数组num和dem则表示数字滤波器系统函数的分子多项式和分母多项式。9) 脉冲响应不变法 numd,dend=impinvar(num,den,Fs)该函数可实现用脉冲响应不变法将模拟滤波器转换为数字滤波器，num和den分别表示模拟滤波器系数函数的分子多项式和分母多项式，Fs是脉冲响应不变法中的抽样频率，单位是Hz。输出变量numd和dend分别表示数字滤波器的系数函数的分子多项式和分母多项式。10) 双线性变换法 numd,dend=bilinear(num,den,Fs)该函数可实现用双线性变换法将模拟滤波器转换为数字滤波器，num和den分别表示模拟滤波器系数函数的分子多项式和分母多项式，Fs是脉冲响应不变法中的抽样频率，单位是Hz。输出变量numd和dend分别表示数字滤波器的系数函数的分子多项式和分母多项式。三、实验内容1. 已知模拟滤波器的传输函数 1) 用Matlab画出模拟滤波器的幅频响应2) 用脉冲响应不变法将转换成数字滤波器的系统函数，取采样时间间隔，画出滤波器的幅频响应。3) 用脉冲响应不变法将转换成数字滤波器的系统函数，取采样时间间隔。画出滤波器的幅频响应。2. 利用Butterworth低通滤波器及脉冲响应不变法设计满足下列指标的数字滤波器：1) 取采样时间间隔，画出滤波器的幅频响应。2) 取采样时间间隔，画出滤波器的幅频响应。3. 用双线性变换法设计一个Chebyshev型数字带通滤波器，要求在100200Hz通带内纹波不大于3dB，通带两边各50Hz外是阻带，衰减不小于40dB，抽样频率为1000Hz。4. 信号，试确定设计指标，设计IIR数字高通滤波器实现滤除的成分。要求利用buttord函数求解滤波器的阶次；利用butter函数设计各IIR数字滤波器；画出滤波器的幅频响应和相频响应；给出IIR数字滤波器的系统函数。四、思考题1. 在实验内容第一题，当采样时间间隔和时，比较两种采样时间下的幅频响应，并分析产生的原因。2. 在实验内容第二题，采用不同的采样间隔，的幅频响应是否一样，为什么，解释其中的原因。3. 双线性变换是一种非线性变换，在实验中你观察到这种现象了吗？应该怎样从哪种数字滤波器的幅频响应中可以观察到这种非线性关系？试举例子说明。实验五 用窗函数法设计FIR数字滤波器实验五 用窗函数法设计FIR数字滤波器一、实验目的：1. 掌握用窗函数法设计FIR数字滤波器的原理和方法2. 熟悉线性相位FIR数字滤波器特性。3. 了解各种窗函数对滤波特性的影响。二、实验原理线性相位特点在实际应用中非常重要，如在数据通信、图像处理、语音信号处理等领域，往往要求系统具有线性相位特性，因而常采用容易设计成线性相位的有限冲激响应FIR数字滤波器来实现。1. 广义线性相位一般来说，一个系统函数可以表示成: （5.1）其中称为系统函数的幅度响应，是一个非负的实函数，而则称为系统函数的相位响应。如果相位响应满足，即系统的群延迟是一个与无关的常数，称系统具有严格线性相位。由于严格线性相位在数学上处理较为困难，在FIR滤波器设计中一般使用广义线性相位。 （5.2）其中和是无关的常数，是一可正可负的实函数，则称系统是广义线性相位，被称为系统的幅度函数2. 线性相位系统频域特性由于的对称性及滤波器阶数的奇偶不同，四种不同类型的线性相位FIR系统的频域特性都有各自的特点。1) 型线性相位滤波器（偶对称，为偶数）型线性相位滤波器的频率响应可以写为 （5.3）系统的幅度函数为 （5.4）型线性相位滤波器的幅度函数的周期为，且关于和偶对称，即 （5.5）2) 型线性相位滤波器（偶对称，为奇数）型线性相位滤波器的频率响应可以写为 （5.6）系统的幅度函数为 （5.7）型线性相位滤波器的幅度函数的周期为，关于偶对称，关于奇对称，且，即 （5.8）所以型线性相位滤波器不能用于高通和带阻等滤波器的设计。3) 型线性相位滤波器（奇对称，为偶数）型线性相位滤波器的频率响应可以写为 （5.9）系统的幅度函数为 （5.10）型线性相位滤波器的幅度函数的周期为，关于和奇对称，且，即 （5.11）所以型线性相位滤波器不能用于高通和低通滤波器的设计。4) 型线性相位滤波器（奇对称，为奇数）型线性相位滤波器的频率响应可以写为 （5.12）系统的幅度函数为 （5.13）型线性相位滤波器的幅度函数的周期为，且关于奇对称，关于偶对称，且即 （5.14）所以型线性相位滤波器不能用于低通滤波器的设计。3. 泄漏现象与Gibbs（吉伯斯）现象以理想滤波器加窗获取FIR滤波器为例来说明泄漏现象与Gibbs（吉伯斯）现象，理想滤波器的截断过程可以表示为： （5.15）其中为所加的窗函数，由DTFT的性质可以得到的频率响应为： （5.16）其中是理想滤波器的频率响应，是窗函数的频率响应，因此，逼近的质量，取决于窗函数的频谱。假设是长度为的矩形窗，则有 （5.17）矩形窗的频谱函数的幅度谱，如图5.1所示 图5.1主瓣处有一个峰值，表示窗函数主要是由直流分量组成。由于矩形窗函数在两个端点的突然截断，使得频谱中存在许多高频分量，这种现象称为频率泄露。如果将理想滤波器的频率响应用幅度函数和相位表示为： （5.18）则可得滤波器的频率响应： （5.19）所以滤波器的幅度函数 （5.20）由上式可以看出来，加窗的滤波器的幅度函数是理想滤波器的幅度函数与窗函数的卷积积分，由于窗函数旁瓣的影响，会形成FIR滤波器的过渡带。这种因为将理想滤波器的单位脉冲响应截断获得的FIR滤波器的幅度响应在通带和阻带都呈现出振荡现象，这个现象就称为Gibbs现象。4. 常用窗函数： 1) 矩形窗 （5.21）2) Hann（汉纳）窗 （5.22）3) Hamming（汉明）窗 （5.23）4) Blackman（布莱克曼）窗 （5.24）5) Kaiser（凯泽）窗 （5.25）其中下面介绍用窗函数设计FIR滤波器的步骤：a) 根据技术要求确定待求滤波器的单位取样响应。b) 根据对过渡带和阻带衰减的要求，选择窗函数的形式，并估计窗函数长度，决定于窗口的形式， 表示滤波器的过渡带。c) 利用选好的窗函数计算滤波器的单位取样