数字滤波器及MATLAB实现.doc

学 号 毕业设计（论文）题目：数字滤波器及MATLAB实现 作 者 届 别 系 别 专 业 指导教师 职 称 完成时间 I摘 要随着信息时代和数字世界的到来，数字信号处理已成为当今一门极其重要的科学和技术。许多的数字信号处理都是把采样频率看成是固定值，但在有的数字系统中，要求数字系统工作在“多采样率”的状态，这就需要改变采样频率，即进行频率转换。本设计研究了变采样率数字滤波器的内插和抽取的原理，并在此基础上，用MATLAB编程设计实现了该变采样率的数字滤波器。在频率转换的设计过程中主要讨论了两个问题，其一是信号的整数倍内插及MATLAB实现；其二是信号的整数倍抽取及MATLAB实现。滤波器的整个设计过程中都是按照理论分析、编程设计、MATLAB实现的步骤来完成。关键词：变采样率；FIR数字滤波器；内插；抽取 ABSTRACTWith the coming of message time and digital world ,the digital signal processing has become one of important science and technology .In most of digital signal processing,the sampling frequency is consider as constant .However multi-sampling frequency is need in some digital system ,thus changing the sampling frequency is unavoidable .namely, frequency convert.In this paper ,the principle of interpolate and decimate of frequency convert digital filter are investigated,and the frequency convert digital filter is realized by programming using MATLAB software. Two major points are discussed about the design of sampling frequency change :the interpolation of the signal of integer and reslization of MATLAB ;the other is about the design of the decimation of the signal of interger and reslization of MATLAB.The design of filter follows the procedures of theoretical analysis ,programming design and concrete realization.Key words: multi-sampling frequency;FIR digital filter;interpolate;decimate 20目 录摘要 ABSTRACT 1 绪论 11.1 数字滤波器的发展概括 11.2 变采样率数字滤波器的研究动态 21.3 研究任务 32 FIR数字滤波器常用设计方法 32.1 MATLAB工具的简介 32.2 FIR数字滤波器线性相位条件 42.3 用窗函数法设计FIR数字滤波器 52.4 频率采样法设计FIR数字滤波器 73 变采样率数字滤波器的基本原理 83.1 信号的整数倍内插 83.2 信号的整数倍抽取 103.3 内插和抽取的综合124 FIR数字滤波器的变采样率实现 124.1 内插的MATLAB实现 134.2 抽取的MATLAB实现 144.3 实际应用举例 154.4 MATLAB实现程序 165 总结 19参考文献 19致谢 201 绪论1.1 数字滤波器的发展概括数字滤波器是数字信号处理理论的一部分。数字信号处理科学中的一项重大的进展是关于数字滤波器设计方法的研究。关于数字滤波器，早在40年代末期，就有人讨论过它的可能性问题，在50年代也有人在研究生班讨论过数字滤波的问题。但直到60年代中期，才开始形成关于数字滤波的一整套完整的正规理论。在这一时期，提出了各种各样的数字滤波器结构，有的以运算误差小为特点，有的则以运算速度高见长，而有的则二者兼而有之；出现了数字滤波器的各种逼近方法，对递归和非递归两类滤波器作了全面的比较，统一了数字滤波器的基本概念和理论。数字滤波器领域的一个重要的进展是对有限冲激响应(FIR)和无限冲激响应(IIR)关系的认识的转化。在初期，一般认为IIR滤波器比FIR滤波器具有更高的运算效率，因而明显地倾向前者，但当人们提出了快速傅立叶变换（FFT）实现卷积运算的概念之后，发现高阶FIR也可以用很高的运算效率来实现，这就促使人们对高性能FIR滤波器的设计方法和数字滤波器的频域设计方法进行了大量的研究，从而出现了这些均属数字滤波器的早期研究。早期的数字滤波器尽管在语音、声纳、地震和医学的信号处理中曾经发挥过作用，但由于当时计算机主机的价格很昂贵严重阻碍了专用数字滤波器的发展。70年代科学技术蓬勃发展，数字信号处理开始与大规模和超大规模集成电路技术、微处理器技术、高速数字算术单元、双极型高密度半导体存储器、电荷转移器件等新技术、新工艺结合起来，并且引进了计算机辅助设计方法，它使数字滤波器的设计仅仅是对相应模拟滤波器的逼近。一般说来，通过对模拟滤波器函数的变换来设计数字滤波器，很难达到逼近任意频率响应或冲激响应，而采用计算机铺助设计则有可能实现频域或时域的最佳逼近，或频域时域的最佳逼近。这样，数字滤波器的分析和设计其内容也更加丰富起来：各种新的数字信号处理系统也都能用专用数字硬件实时加以实现。所谓数字滤波器，是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些成分的器件。因此，数字滤波的概念和模拟滤波相同，只是信号的形式和实现滤波方法不同。正因为有该不同点，数字滤波器具有比模拟滤波器的精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及实现模拟滤波器无法实现的特殊滤波功能。滤波器的种类很多，从功能上可分为低通、高通、带通和带阻滤波器，每一种又有模拟滤波器和数字滤波器两种形式。如果滤波器的输入和输出都是离散时间信号，则该滤波器的冲击响应也必然是离散的，这种滤波器称之为数字滤波器。该滤波器通过对时域中离散的采样数据作差分运算实现滤波。一个线性时不变数字滤波器可以用常系数线性差分方程表示为：(1.1) 式中x(n),y(n)分别是输入和输出信号序列；分别是滤波系数。从实现的网络结构或者从单位脉冲响应分类上可分为无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。当(1)式中系数全部为零时，就有 (1.2)这种形式的滤波器为有限冲击响应滤波器，简称FIR型，此时系统的输出仅与输入有关，它的实现一般采用非递归的算法。1.2 变采样数字滤波器的研究动态随着信息时代和数字世界的到来，数字信号处理已成为当今一门极其重要的学科和技术领域。数字信号处理在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理应用中，数字滤波有广泛的应用。 我们以前接触到的数字信号处理大都是把采样频率看成是固定值的情况，但是在实际系统中，会经常遇到要求数字系统能够工作在“多采样率”的状态，这就要求我们按照需要来改变采样频率，即进行频率转换。一般认为，在满足采样定理的前提下，先将以采样率F1采集的数字信号进行D/A转换成模拟信号，再按采样率F2进行A/D变换，来实现从F1到F2的采样率转换。但是这样很容易使信号受到伤害，所以我们不经过数模转换，直接在数字域中实现。在数字域中，我们通过“抽取”来降低采样率去掉多余的数据，通过“插值”来提高采样率增加我们需要的频率，来实现采样率的转换。频率转换技术也称过采样技术，既可以用于不同速率的信号处理，也可以应用于A/D，D/A转换器，目的是将模拟滤波器的部分指标由数字滤波器承担，使A/D前的抗混叠滤波及D/A后的平滑滤波器容易实现。变采样率数字信号处理在数字滤波设计、序列编码和数据压缩等很多领域得到了广泛的应用，与滤波器组理论和小波变换也有密切联系，现在很多实际应用中都要用到采样率的转换（变采样率），比如，数字电视系统、数字电话系统、对一个非平稳随机（语音信号）作谱分析或编码等等方面。而我所要做的就是利用MATLAB编程实现来进行抽取和内插，从而实现采样率的转换。通常我们在FIR中实现多采样率转换。MATLAB是一套用于科学计算的可视化高性能语言与软件环境。它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个界面友好的用户环境。它的信号处理工具箱包含了各种经典的和现代的数字信号处理技术，是一个非常优秀的算法研究与辅助设计的工具。在设计数字滤波器时，通常采用MATLAB来进行辅助设计和仿真。MATLAB作为一种科学计算软件，为数字信号处理系统提供了很多快捷的计算和仿真工具。MATLAB信号处理箱为我们这个研究课题提供了一个变采样率函数resample,用它可以把变采样率的全过程在一步中完成.即我们通过MATLAB编程来实现采样率的转换。1.3 设计任务本设计对实现变采样率数字滤波器的内插和抽取方法作了研究，并在此基础上用MATLAB编程实现变采样率的数字滤波器。2 FIR数字滤波器常用设计方法2.1 MATLAB工具的简介MATLAB是由美国MathWorks公司开发成功的一种用于科学工程计算的高效率高级语言。MATLAB的语法规则简单，很贴近人的思维方式。用MATLAB编写程序，犹如在一张演算纸上排列公式和求解问题一样效率高，因此被称为“演算纸式的”科学工程算法语言。MATLAB的基本数据单位是矩阵，它的命令表达与数学、工程中常用的习惯表达十分相似。比如，矩阵方程b=Ax,在MATLAB中被写成b=A*x。而若要通过A、b和x,那么只要写出命令x=Ab即可，而完全不需要对矩阵的乘法和求逆进行编程。在MATLAB环境下，许多复杂的数学运算，如求矩阵的行列式值、求函数的微分、求函数的积分、进行多项式插值、解微分方程等，都有现成的库函数可以调用，并且各个库函数是根据不同的应用情况采用不同的优化算法的，保证了结果的可靠性和求解的快速性。经过MathWorks公司的不断完善，MATLAB版本不断更新，其功能也越来越强，使她在诸如一般数值计算、数字信号处理、系统识别、自动控制、振动理论、时序分析与建模、优化设计、神经网络控制、化学统计学、动态仿真系统、特殊函数和图形领域表现出一般高级语言难以比拟的优势，可以方便地用于几乎所有的科学和工程计算的各个方面。可以说，MATLAB不仅是一种编程语言，而且在广义上是一种语言开发系统。MATLAB的语言特点如下：1. 高效方便的矩阵和数组运算；2. 编程效率高；3. 结构化/面向对象；4. 方便的绘图功能；5. 用户方便使用；6. 强大的工具箱；7. 扩充能力强。MATLAB还提供了SIMULINK工具箱，大大方便了复杂控制系统的仿真，便于系统分析。SIMULINK模型的仿真包含一系列一般微分方程的数值积。SIMULINK为仿真提供了许多种计算这些方程的求解器。由于不同动态系统的行为差异，所以在解决特定的问题时，某些求解器会比其它的更有效一些。为得到精确且迅速的计算结果应选择合适的求解器，并且设定合适的参数。求解器可以选择变步长或是定步长。变步长的求解器在仿真期间会改变其步长，它们提供误差控制和过零点检测。定步长的求解器在仿真期间采用相同的步长，它们不提供误差控制也不定位过零点。仿真性能和精度由多种因素决定，包括模型的设计和仿真参数的选择。求解器使用它们的缺省参数值可以使大多数模型的仿真比较精确有效，然而，对于一些模型如果调整求解器和仿真参数将会产生更好的结果，而且，如果对模型的性能比较熟悉，并且将这些信息提供给求解器，得到的仿真效果将会提高。仿真程序的运行可通过命令运行也可通过菜单运行仿真。通过命令行运行仿真与通过菜单运行仿真相比，有许多优点，既可以模仿M文件和MEX文件形式，甚至是SIMULINK模块图形式的模型；也可以运行M文件来运行仿真：这样，仿真和模块参数可以反复地更改。在MATLAB的命令窗口或M文件中，输入仿真命令也可以运行仿真。在进行MonteCarlo分析时，可任意地改变参数，并在一个循环中分别运行仿真，得到不同参数的结果。可在命令行中输入sim和set param命令来运行仿真。2.2 FIR数字滤波器线性相位条件根据数字滤波器冲激响应的时域特征，数字滤波器分为两种，即无限长冲激响应滤波器(IIR)和有限长冲激响应滤波器(FIR)。FIR有如下自己突出的优点：系统总是稳定的、易于实现线性相位、允许设计多通带(或多阻带)滤波器。因此FIR在数字信号处理中得到广泛的应用。但与IIR相比，在满足同样的阻带衰减的情况下需要较高的阶数。滤波器阶数越高将占用更多的DSP运算时间。因此，对FIR的设计目标是在满足指标要求的情况下尽量减少滤波器的阶数。数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。FIR的冲激响应h(n)是有限长的，M阶FIR可以表示如下： (2.1)系统函数为： (2.2)在本设计中，我们采用满足第一类线性相位的FIR滤波器，它的幅度函数为： (2.3)FIR数字滤波器具有线性相位特点的条件是它的单位冲激响应h(n)满足偶对称，即h(n)=h(N-n-1) 或h(n)=-h(N-n-1)其中：NFIR滤波器的长度。若h(n)满足偶对称，则其相频特性为： (2.4)若h(n)满足奇对称，则其相频特性为： (2.5)从频域上看：h(n)是长度为N的实序列，所以h(n)的N点DFT满足共轭对称性，即： (2.6)而幅度函数则一定呈偶对称特性，即： (2.7)它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。在本设计中用窗函数法。2.3 用窗函数法设计FIR数字滤波器设希望设计的滤波器传输函数为，是与其对应的单位脉冲响应 ，因此， (2.8) (2.9) 如果能够由已知的求出，经过Z变换可得到滤波器的系统函数。但一般情况下，逐段恒定，在边界频率处有不连续点，因而是无限时宽的，且是非因果序列，例如，理想低通滤波器的传输函数为：(2.10)相应的单位取样响应为： (2.11)由上式看到，理想低通滤波器的单位取样响应是无限长，且是非因果序列。为了构造一个长度为 N的线性相位滤波器，只有将截取一段，并保证截取的一段对(N-1)/2对称。设截取的一段用h(n)表示，即 (2.12)式中是一个矩形序列，长度为N。当a取值为(N-1)/2时，截取的一段h(n)对(N-1)/2对称，保证所设计的滤波器具有线性相位。我们实际实现的滤波器的单位取值响应为h(n),长度为N,其系统函数为H(z), (2.13)这样我们用一个有限长的序列h(n)去代替，肯定会引起误差，表现在频域就是通常所说的吉布斯效应。该效应引起通带内和阻带内的波动性，尤其使阻带的衰减小，从而满足不了技术上的要求。这种吉布斯效应是由于将直接截断引起的，因此也称截断效应。用窗函数法设计滤波器的程序流程图如图2.1所示。 读入窗口长度N开始计算调用窗函数子程序计算调用子程序（函数）计算调用绘图子程序（函数）绘制幅度相位曲线结束图2.1 用窗函数法设计滤波器程序流图2.4 频率采样法设计FIR数字滤波器窗函数法是以h(n)为媒介的时域设计法，而滤波器指标往往是在频域给出的，为此，要由H算出，加窗后又从h(n)算出来检验。当理想频响H是任意曲线，或者不存在明确的解析表达式时，求就困难一些。因此，我们不得不想：能否不要频域-时域-频域地反复，而直接从频域来设计呢？这就是采用FIR滤波器的频域设计方法频域采样法。频率采样法先对理想频响采样，得到样值H(k),再利用插值公式直接求出系统转换函数H(z)，以便实现；或者求出频响H,以便与理想频响进行比较。在区间上对进行N点采样，等效于时域以N为周期延迟。频率采样法的步骤可归纳为：1. 给定理想频响。2. 确定采样点数，对理想频响采样得到H(k)。3. 代入下式，即得FIR数字滤波器的转移函数： (2.14)频率采样法可以看作为插值法，它在采样的上保证等于而在非插值点(采样点)上，是插值函数的线性叠加。这种方法的缺点是通带和阻带的边缘要精确确定。窗口法的矩形窗是一种最小平方逼近法，它不能保证在每个局部位置误差都很小。3 变采样率数字滤波器的基本原理我们以前接触到的数字信号处理大都是把采样频率看成是固定值的情况，但是在实际系统中，会经常遇到要求数字系统能够工作在“多采样率”的状态，这就要求我们按照需要来改变采样频率，即进行频率转换。采样率转换处理是一种用于改变已经数字化的离散信号有效采样率的一种有用而神奇的方法。采样率转换有很多用处，当信号通过低通滤波以后带宽变窄时，通过降低采样率可以减少计算的工程量。在实时的处理中，当两种分开的硬件处理器按两中不同的采样率运算时，其数字信号的数据必须互相转换，因而不得不进行采样率的转换。在卫星和医学图象处理中，采样率转换对于图象增强，尺度改变以及图象旋转等处理也是必须的，采样率转换还可以用于减少某些窄带数字滤波器的复杂性。一般认为，在满足采样定理的前提下，先将以采样率F1采集的数字信号进行D/A转换成模拟信号，再按采样率F2进行A/D变换，来实现从F1到F2的采样率转换。但是这样很容易使信号受到伤害，所以我们不经过数模转换，直接在数字域中实现。多采样率数字信号处理是指在一个系统中有2个或2个以上的采样速率，通过有目的地改变采样速率，以便于信号的处理，达到节省计算量的目的。在数字域中，我们通过“抽取”来降低采样率去掉多余的数据，通过“插值”来提高采样率增加我们需要的频率，来实现采样率的转换。抽取是降低采样率以去掉多余数据的过程，而插值则是提高采样率已增加数据的过程。3.1 信号的整数倍内插整数倍内插是在已知的相邻两个采样点之间插入I-1个样本点。由于这些样本点并非已知的值，所以关键问题是如何求出着I-1个样本值。前提是不经过模拟域，要直接在数字域完成这个插值，也称之为升采样过程。用数字方法实现整数内插的步骤如下：先在已知采样序列x()的相邻两个样点之间等间隔插入I-1个零值样本点，假如输出为y，则输入输出的关系为：(3.1)式中I为大于1的整数，称为插入因子。然后经过低通滤波，把滤波得到的结果带入这些位置，得到内插的结果。信号内插后，序列的频谱在频率方向压缩了I倍。由于原数字频谱假设为的周期函数，故压缩 之后的频谱将在0范围内产生I-1个镜像频谱。内插结果虽然没有造成频谱混叠，但频谱压缩之后形成的 I-1个镜像频谱并没有带来任何有效信息，必须清除。因此，必须在内插之后增加一级数字低通抗镜像滤波。零值的内插如图3.1所示。图3.1 零值内插方案的系统框图其中 ， 为：原信号的采样序列； 为：内插之后的采样序列； 为：低通滤波器； 为：对内插后采样序列进行滤波得到的序列。内插系统如图3.2所示。图3.2 数字信号的整数倍内插MATLAB程序实现如下：nx=0：Nx-1; %设定x的自变量向量x=exp（-alpha*nx*T）; %给出x序列的值Ny=Nx*I;ny=0：Ny-1； %设定y的自变量向量P=ones（1，I*Nx）； %生成频率为I*Nx的单位脉冲序列py=zeros（1，I*Nx）; % 先y序列初始化（全置零）y（1：I：I*Nx）=x b=fir1(N,wc) %滤波器系数y1=filter（b,1,y）; %用滤波器进行滤波y2=exp（-alpha*ny*T/I）/I； %原始函数用高采样率的序列X=fftshift（fft（x）; %求x频谱，并移到对称位置Y=fftshift（fft（y）; %求y频谱，并移到对称位置p=fftshift（fft（p）; %求p频谱，并移到对称位置Y1=fftshift（fft（y1）; %求y1频谱，并移到对称位置Y2=fftshift（fft（y2）; %求y2频谱，并移到对称位置3.2 信号的整数倍抽取整数倍抽取是在已知序列x中每隔D-1个采样点取出一个样本点，组成新的序列y的过程。因此y序列的所有值都是已知的，可以直接写出下列表达式： (3.2)也就是说，输入序列中下标为D的整数倍的点上的x都被取入y序列，而其他的x序列中的样本都被舍去了。信号抽取后，序列的频谱在频率方向扩张了D倍。由于数字频谱假设是的周期函数，故扩展后的频谱将在0范围内可能产生频谱混叠。事实上，若原序列采样满足采样定理，抽取采样将采样频率将至，就有可能破坏采样条件，从而导致频谱混叠。为原信号的采样频率。为了避免抽取采样之后频谱在0范围内产生混叠，要求原序列的数字频谱控制在0之间。这样，抽取采样之后的最大扩展带宽为0，不会造成混叠。为了达到这个目的，在抽取采样之前一般需要增加一级数字低通抗混叠滤波。抽取方案如图3.3所示。图3.3 抽取方案的系统框图其中 ， 为：原信号的采样序列； 为：低通滤波器； 为：对原信号进行滤波得到采样序列； 为：对经过滤波的信号进行抽取得到的序列。抽取系统如图3.4所示。图3.4 数字信号的整数倍抽取MATLAB程序实现如下：nx=0:Nx-1;x=exp(-alpha*nx*T); %给出x序列的值b=fir1(N,wc) %滤波器系数x1=filter(b,1,x);Ny=ceil(Nx/M);ny=0:Ny-1;Y1=x1(1:D:Nx); %给出y1序列在x1对应点处的值 X=fftshift（fft（x）; %求x频谱，并移到对称位置X1=fftshift（fft（x1）; %求y频谱，并移到对称位置Y1=fftshift（fft（y1）; %求y1频谱，并移到对称位置 3.3 抽取和内插的综合虽然按整数倍因子抽取或按整数因子内插改变采样率很有用，但如果我们需要按非整数因子改变采样率怎么办？现在有一个好办法，可以按有理数I/D实现采样率转换，即先按整数因子I插值，再按整数因子D抽取。因为比值I/D可以按我们希望的精度得到，所以通过正确选择整数I和D，我们实际上可按几乎任何因子改变采样率。本文主要介绍整数倍的内插和抽取。4 FIR数字滤波器的变采样率实现由于FIR结构的绝对稳定性且很容易做成线性相位，特别是容易实现高效结构，所以用FIR结构实现多采样率系统具有很大的优越性，故在多采样率系统的实现中绝大多数采用FIR结构。程序执行框图如4.1所示。 图4.1 程序的执行框图输入判断语句k=0对信号进行内插抗镜像低通滤波绘制时域和频域波形结束开始Y抗混叠滤波N对信号进行抽取绘制时域和频域的波形4.1 内插的MATLAB实现当I=2时的内插如图4.2所示。图4.2 当I=2时信号的内插当I=3时的内插如图4.3所示。图4.3 当I=3时信号的内插 4.2 抽取的MATLAB实现当D=3时的抽取如图4.4所示。图4.4 当D=3时信号的抽取当D=2时的抽取如图4.5所示。图4.5 当D=2时信号的抽取4.3 实际应用举例我们将一个给定信号x=exp(-5*nx/fs)，将其采样频率由4kHZ变为12kHZ,实现如图4.6所示。图4.6 信号的升采样频率我们将一个给定信号x=exp(-5*ny/fs)，而将其采样频率由8kHZ变为4kHZ,实现如图4.7所示。图4.7 信号的降采样频率4.4 MATLAB实现程序MATLAB的实现程序为：N=input;n=0:N;wc=input;alpha=(N-1)/2;m=n-alpha+eps;hd=sin(wc*m)/(pi*m);w_ham=(hamming(N);h=hd.*w_ham;theta=-alpha*wc;H=h*cos(n-(N-1)/2)*wc);k=1;if k=0;nx=0:N-1;T=input; x=exp(-alpha*nx*T);I=input;Ny=N*I;ny=0:(Ny-1);y=zeros(1,I*N);y(1:I:I*N)=x;b=fir1(N,wc);y2=filter(b,1,y);X=fftshift(fft(x);Y=fftshift(fft(y);Z=fftshift(fft(y2);X1=fft(x);Y1=fft(y);Z1=fft(y2);figure(3)subplot(3,2,1);stem(abs(X);title(原信号序列);k1=0:1:nx;w1=2*pi/(nx+1)*k1;subplot(3,2,2);plot(w1/pi,abs(X1);title(原信号频谱);subplot(3,2,3);stem(abs(Y);title(内插序列);k2=1:1:N;w2=4000/(N+1)*k2;subplot(3,2,4);plot( w2/pi,abs(Y1);title(内插频谱)subplot(3,2,5);stem(abs(Z);title(滤波后序列);k3=0:1:ny;w3=4000/(ny+1)*k3;subplot(3,2,6);plot(w3/pi,abs(Z1);title(滤波后频谱);elsefigure(3)nx=0:N-1;T=input;x=exp(-alpha*nx*T);b=fir1(N,0.3*pi);x1=filter(b,1,x)D=input;Ny=ceil(N/D);ny=0:Ny-1;y=x1(1:D:N);y2=exp(-alpha*ny*T*D)*D;z=fft(y2);X=fftshift(fft(x);X1=fft(x);X1=fftshift(fft(x1);Y=fftshift(fft(y);Y1=fft(y);Z=fftshift(fft(y2);Z1=fft(y2);subplot(3,2,1);stem(abs(X);title(原信号序列);k1=0:1:nx;w1=2*pi/(nx+1)*k1;subplot(3,2,2);plot(w1/pi,abs(X1);title(原信号频响);subplot(3,2,3);stem(abs(Y);title(抽取序列);k2=1:1:N;w2=2*pi/(N+1)*k2;subplot(3,2,4);plot(w2/pi,abs(Y1);title(抽取频响);subplot(3,2,5);stem(a