基于MATLAB函数直接实现椭圆滤波器的理论设计及应用分析样本.doc

本科课程论文基于MATLAB函数直接实现椭圆滤波器理论设计及应用分析 专 业： 指引教师： 学生姓名： 学生学号： 中华人民共和国贵州贵阳6月目录摘要III引 言4第一章 椭圆滤波器基本理论51.1椭圆滤波器概述51.2椭圆滤波器设计数学推导51.3关于归一化讨论8第二章 椭圆滤波器设计92.1椭圆滤波器设计构造图92.2设计椭圆数字滤波器环节92.3数字椭圆低通滤波器MATLAB实现92.3.1设计椭圆滤波器所用函数92.3.2频谱分析所用函数102.3.3 其她命令函数122.4 椭圆低通滤波器设计程序12第3章 仿真图133.1原始信号及其采样仿真图，如图4所示133.2信号通过椭圆低通滤波器仿真图，如图5所示133.3 椭圆低通滤波器幅度特性，如图6所示143.4 对信号进行傅里叶变换仿真，如图7所示14第四章 椭圆滤波器在语音去噪中应用分析154.1 语音信号采集154.2语音信号频谱分析164.3 设计椭圆滤波器184.4 信号滤波解决19心得体会22参照文献23附录24附录一：椭圆低通滤波器设计程序：24附录二：第四章源程序清单25基于MATLAB函数直接实现椭圆滤波器理论设计及应用分析摘要近代电信装备和各类控制系统中，滤波器应用极为广泛；在所有电子部件中，使用最多，技术最复杂要算滤波器了。滤波器优劣直接决定产品优劣，因此对滤波器研究和生产从来为各国所注重。随着当代科学技术发展，滤波器在咱们研究中占着越来越大份额，它影响真咱们信号技术研究与发展，滤波器所带来巨大影响和作用使咱们有必要去探讨它应用和发展。滤波器是一种能使有用信号顺利通过而同步对无用频率信号进行抑制（或衰减）电子装置。工程上惯用它来做信号解决、数据传送和抑制干扰等。滤波器模仿滤波器和数字滤波器之分。模仿滤波器由有源和无源之分，无源滤波器重要是R，L，C构成。模仿滤波器会有电压漂移，温度漂移和噪声等问题。搭建模仿滤波器和数字滤波器之间桥梁是采样定理，采样定理将持续信号转化成数字信号。模仿滤波器特性可以用其频率响应来描述，按其特性不同，可以分为低通滤波器，高通滤波器，带通滤波器和带阻滤波器等。本文将通过运用MATLAB滤波滤波器设计函数直接实现椭圆滤波器设计，找到应用MATLAB来设计椭圆滤波器办法。简介了椭圆型滤波器基本理论和设计思想，给出了基于MATLAB设计低通、带通、高通椭圆型滤波器详细环节和运用MATLAB产生一种包括低频、中频、高频分量持续信号，并实现对信号进行采样。文中还对采样信号进行频谱分析和运用设计椭圆滤波器对采样信号进行滤波解决，并对仿真成果进行分析和解决。详细简介了在基于MATLAB设计椭圆滤波器过程中惯用到工具和命令。核心字：低通，滤波器，MATLAB，持续信号引 言 信号解决是科学研究和工程技术许多领域都需要进行一种重要环节，老式上对信号解决大都采用模仿系统实现。随着人们对信号解决规定日益提高，以及模仿信号解决中某些不可克服缺陷，对信号许多解决而采用数字办法进行。近年来由于大规模集成电路和计算机技术进步，信号数字解决技术得到了飞速发展。数字信号解决系统无论在性能、可靠性、体积、耗电量、成本等诸多方面都比模仿信号解决系统优越多，使得许多以往采用模仿信号解决系统越来越多地被数字解决系统所代替，数字信号解决技术在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛应用。在数字信号解决中，数字滤波器十分重要并已获得广泛应用，数字滤波器与模仿滤波器比较，具备精度高、稳定、体积小、重量轻、灵活、不规定阻抗匹配以及实现模仿滤波器无法实现特殊滤波功能等长处。在各种滤波器中，椭圆滤波器具备其独特长处。本次设计中所用到数学软件为MATLAB。MATLAB和Mathematica、Maple并称为三大数学软件，它是由美国mathworks公司发布重要面对科学计算、可视化以及交互式程序设计高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统建模和仿真等诸多强大功能集成在一种易于使用视窗环境中，为科学研究、工程设计以及必要进行有效数值计算众多科学领域提供了一种全面解决方案，并在很大限度上挣脱了老式非交互式程序设计语言（如C、Fortran）编辑模式，代表了当今国际科学计算软件先进水平。第一章 椭圆滤波器基本理论1.1椭圆滤波器概述惯用数字滤波器类型有巴特沃斯(Butterworth) ，切比雪夫(Chebyshev) 及椭圆型 滤波器，其中椭圆滤波器（Elliptic filter）又称考尔滤波器（Cauer filter），是一种性能优越滤波器。从传递函数来看，巴特沃斯和切比雪夫滤波器传播函数都是一种常数除以一种多项式， 为全极点网络， 仅在无限大阻带处衰减为无限大， 而椭圆函数滤波器在有限频率上既有零点又有极点。极零点在通带内产生等波纹， 阻带内有限传播零点减少了过渡区， 可获得极为陡峭衰减曲线。也就是说在阶数相似条件下，椭圆滤波器相比于其她类型滤波器，能获得更窄过渡带宽和较小阻带波动， 就这点而言， 椭圆滤波器是最优。它陡峭过渡带特性是用通带和阻带起伏为代价来换取，并且在通带和阻带波动相似，这一点区别于在通带和阻带都平坦巴特沃斯滤波器，以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹切比雪夫滤波器。 总结起来，椭圆滤波器具备如下特点：1）椭圆低通滤波器是一种零、极点型滤波器，它在有限频率范畴内存在传播零点和极点。2）椭圆低通滤波器通带和阻带都具备等波纹特性，因而通带，阻带逼近特性良好。 3）对于同样性能规定，它比前两种滤波器所需用阶数都低，并且它过渡带比较窄。但是椭圆滤波器传播函数是一种较复杂逼近函数， 运用老式设计办法进行电路网络综合要进行繁琐计算， 还要依照计算成果进行查表， 整个设计， 调节都十分困难和繁琐。而用MATLAB设计椭圆滤波器可以大大简化设计过程。1.2椭圆滤波器设计数学推导椭圆滤波器振幅平方函数为 ： （1）其中是雅可比(Jacobi) 椭圆函数，雅可比椭圆函数是阶数N有理函数，N=5时特性曲线如图1所示。图1 N=5时雅可比椭圆函数特性曲线 由图1 可见，在归一化通带内( - 1 1) ，() 在(0 ，1) 间振荡，而超过 后，() 在( ， ) 间振荡。L 越大，也变大。这一特点使滤波器同步在通带和阻带具备任意衰减量。L 是一种表达波纹性质参量。雅可比椭圆函数还具备如下性质： （2）阶数N等于通带和阻带内最大点和最小点总和，为与通带衰减关于参数。系统函数和阶数N是由系统下面性能指标来拟定，重要有：截止频率，通带内最大衰减和阻带截止频率以及阻带内最小衰减。假定是频率归一化基准频率，即 （3）定义频率选取性因数k为 （4）则截止频率分别归一化为 （5）再次假定 （6） （7） （8） （9）则得到椭圆滤波器阶数N为 （10）这时，令归一化基准频率为，则得到归一化后椭圆低通滤波器系统函数为 （11）式中， 因此，实际椭圆低通滤波器就可以由归一化系统函数来得到 （12） 图2 为典型N 为奇数椭圆滤波器幅度特性， 当，和A 拟定后，阶次N 即可拟定，进而可以设计出椭圆滤波器。图2 椭圆滤波器幅度特性1.3关于归一化讨论归一化是一种简化计算方式，重要是为了数据解决以便提出来，即将有量纲表达式，通过变换，化为无量纲表达式，成为纯量。例如，复数阻抗可以归一化书写：Z = R + jL = R(1 + jL/R) ，复数某些变成了纯数量了，没有量纲。归一化办法（Normalization Method）把数据映射到01范畴之内解决，更加便捷迅速，应当归到数字信号解决范畴之内。其详细作用是归纳统同样本记录分布性。归一化在01之间是记录概率分布，归一化在-1+1之间是记录坐标分布。归一化有同一、统一和合一意思。无论是为了建模还是为了计算，其基本度量单位要同一。但是归一化解决并不总是适当，依照输出值分布状况，原则化等其他记录变换办法有时也许更好。详细状况还需详细分析。第二章 椭圆滤波器设计2.1椭圆滤波器设计构造图椭圆滤波器设计构造图如图3所示：图3 构造框图2.2设计椭圆数字滤波器环节由于模仿滤波器设计办法非常成熟， 许多典型系统有成熟公式、图表可以查阅，便于设计；因而设计数字滤波器重要办法是：一方面设计一种适当模仿滤波器， 然后将她“ 变换”成满足给定指标数字滤波器。设计椭圆数字滤波器环节:(1) 拟定数字滤波器性能指标、;(2) 将数字滤波器性能指标转换成相应模仿滤波器性能指标;(3) 设计满足指标规定模仿滤波器;(4) 通过变换将模仿滤波器转换成数字滤波器。2.3数字椭圆低通滤波器MATLAB实现Matlab 是MathWorks 公司于1984 年正式推出一套集数值计算、符号运算及图形解决等强大功能于一体科学计算语言。作为强大科学计算平台，她几乎可以满足所有计算需求。她应用范畴涵盖了当今几乎所有领域，如电子、半导体制造、医学研究、航空航天、汽车制造、分子模型、影视制作、建筑等行业。Matlab 具备如下优势和特点:和谐工作平台和编程环境，简朴易用程序语言，强大科学计算及数据解决能力，出众图形解决功能，应用广泛模块集和工具箱，实用程序接口和发布平台，模块化设计和系统级仿真。随着Matlab不断完善， 特别是Matlab 信号解决工具箱( SignalProcessing Toolbox) 推出，如今Matlab 已经成为数字信号解决DSP(Digital Signal Processing) 应用中分析和仿真设计重要工具。2.3.1设计椭圆滤波器所用函数Matlab 信号解决工具箱提供了设计椭圆滤波器函数:ellipord 函数、ellip 函数和ellipap函数。1. ellipord 函数功能是求滤波器最小阶数和截止频率，其调用格式：N， = ellipord( ， ， ， ) 可以得到数字椭圆型滤波器最小阶数N和截止频率 ，并使滤波器在通带内(0 ，) 波纹系数不大于通带最大衰减 ，阻带内( ，1) 波纹系数不不大于阻带最小衰减。其中是椭圆滤波器通带截止角频率，是椭圆滤波器阻带起始角频率。依照本次任务书设计规定，需要产生一种持续信号，包括低频5Hz，中频15Hz，高频30Hz三个分量，并对其进行采样，采样频率为100Hz，采样点数为100。设计低通滤波器对信号进行滤波解决，滤除中频和高频信号。由于已知参数有限，对于设计中所用到参数可取= 0.1，=40，通带截止频率Wp=5Hz，阻带截止频率 Ws=10Hz，归一化解决wp=2*Wp/Fs； ws=2*Ws/Fs 。依照程序：Wp=5;Ws=10;Fs=100;rp=0.1;rs=40;wp=2*Wp/Fs;ws=2*Ws/Fs;n，wn=ellipord(wp，ws，rp，rs)可得出：n =4wn =0.10002. ellip 函数功能是设计滤波器，其调用格式:b，a = ellip ( N， ， ，) ， 运用ellipord 函数得到最小阶数N和截止频率，可以设计低通滤波器。其中，b、a分别为椭圆滤波器传播函数分子、分母多项式。3. ellipap函数功能是直接返回椭圆滤波器零点z、极点p和增益k，其调用格式：z ， p ， k=ellipap(N ， )2.3.2频谱分析所用函数Matlab信号解决工具箱提供了频谱分析函数:fft函数、filter函数和freqz函数。1. fft函数功能是对信号进行迅速傅里叶变换，其调用格式： Y = fft(X) Y = fft(X，n)Y = fft(X，dim) Y = fft(X，n，dim)matlabfft序号是从1到n，大多数采用从0到n-1，Y=fft（x）之后，这个Y是一种复数，它模值应当除以（length(x)2），才干得到各个频率信号实际幅值。fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴。由此可以懂得FFT变换数据对称性。因而用FFT对信号做谱分析，只需考察0Nyquist频率范为内福频特性。若没有给出采样频率和采样间隔，则分析普通对归一化频率01进行。此外，振幅大小与所用采样点数关于，例如：采用128点和1024点相似频率振幅是有不同体现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致。为了与真实振幅相应，需要将变换后成果乘以2除以N。2. Freqz函数功能是用来求幅频响应，其调用格式：h，w=freqz(b，a，n)h，f=freqz(b，a，n，Fs)h=freqz(b，a，w)h=freqz(b，a，f，Fs)freqz(b，a，n)阐明： freqz 用于计算数字滤波器H(Z)频率响应函数H(ej)。h，w=freqz(b，a，n)可得到数字滤波器n点复频响应值，这n个点均匀地分布在0，上，并将这n个频点频率记录在w中，相应频响值记录在h中。规定n为不不大于零整数，最佳为2整多次幂，以便采用FFT计算，提高速度。缺省时n =512。h，f=freqz(b，a，n，)用于对在0，/2上等间隔采样n点，采样点频率及相应频响值分别记录在f 和h中。由顾客指定（以HZ为单位）值。h=freqz(b，a，w)用于对在0，上进行采样，采样频率点由矢量w指定。h=freqz(b，a，f，Fs) 用于对在0，上采样，采样频率点由矢量f指定。freqz(b，a，n) 用于在当前图形窗口中绘制幅频和相频特性曲线。3. filter函数功能是运用IIR滤波器和FIR滤波器对数据进行滤波，其调用格式： y，zf=filter(b，a，x) y=filter(b，a，x，zi) y=filter(b，a，x)阐明：filter采用数字滤波器对数据进行滤波，其实现采用移位直接型构造，因而合用于IIR和FIR滤波器。滤波器系统函数为 即滤波器系数a=.，b=.，输入序列矢量为x。这里，原则形式为=1，如果输入矢量a时，1，则MATLAB将自动进行归一化系数操作；如果=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)除得到矢量y外，还得到x最后状态矢量zf。2.3.3 其她命令函数在设计过程中出以上功能函数外，还用到了诸多其她函数。例如： clc 清除命令窗口中内容 Clear 清除内存中变量和函数 Plot 绘制线性图形 Abs 取模 等。2.4 椭圆低通滤波器设计程序参见附录第3章 仿真图3.1原始信号及其采样仿真图，如图4所示图4 原始输入信号及其采样图3.2信号通过椭圆低通滤波器仿真图，如图5所示图5 信号通过椭圆低通滤波器仿真图3.3 椭圆低通滤波器幅度特性，如图6所示图6 椭圆低通滤波器幅度特性3.4 对信号进行傅里叶变换仿真，如图7所示图7 信号傅里叶变换第四章 椭圆滤波器在语音去噪中应用分析语音信号是一种非平稳信号,人们在语音通讯过程中会受到来自周边环境,和传播介质影响,产生噪音,影响人们听觉,因而咱们需要对语音信号进行去噪解决。使用数字滤波器可以有效地去除语音信号中高频和低频噪声,本文重要研究椭圆滤波器在语音去噪中应用。随着科学技术发展,人和机器交流成为各国研究新课题,其中语音辨认是最为重要一种环节。机器在接受语音时候,往往会受到环境噪声和其她噪声影响,因而语音去噪是语音辨认前提。一方面,咱们先理解一下语音信号,语音信号频谱覆盖在50Hz4kHz,较为丰富信号重要集中在1kHz附近,语音信号具备短时平稳性,因而语音信号常被分段或分帧来解决,普通每秒帧数约为33100。因此普通滤波器去噪时必要考虑语音信号自身特性。当前在国内外在语音去噪方面提出了许多去噪办法:自有关相减法、谐波增强法、自适应噪声滤波法等等。本文重要是基于椭圆滤波器,针对详细语音信号研究去噪办法。滤波器是一种选频装置。它对某一种频率或几种频率范畴内信号给以极小衰减,使这某些信号可以顺利通过;对其她频带内信号则给以很大衰减,从而尽量制止这某些信号通过。在各种滤波器中,椭圆滤波器具备其独特长处。4.1 语音信号采集运用Windows下录音机，录制语音信号“人们好，我是XX”，时间在3s左右。然后在Matlab软件平台下，运用函数wavread对语音信号进行采样，记住采样频率和采样点数。录音机使用如图4.1所示图4.1 录音机使用图4.2语音信号频谱分析一方面画出语音信号时域波形；然后对语音号进行迅速傅里叶变换，得到信号频谱特性。 y,fs,bits = wavread(cc.wav); sound (y,fs,bits); plot(y);title(时域波形) t=(1:16000)/8000; plot(t,y);xlabel(t);ylabel(y);title(时域波形) x=x；y=x+sin(fn*2*pi*t)； sound(y,fs,bits)； % 应当可以明显听出有尖锐单频啸叫声其中，y,fs,bits=wavread(cc.wav)； % 输入参数为文献全途径和文献名，输出第一种参数是每个样本值，fs是生成该波形文献时采样率，bite是波形文献每样本编码位数。sound(y,fs,bits)； % 按指定采样率和每样本编码位数回放由绘图命令可以得届时域波形，如图4.2-1所示图4.2-1 时域分析 Y=fft(y); magY=abs(Y);angY=angle(Y); w=(1:16000)/16000*2*pi; plot(w/pi,magY)由绘图命令可得频域幅度谱，如图4.2-2所示图4.2-2 频域幅度谱 f=w/(2*pi)*Fs; f=f(1:8000);magY=magY(1:8000);angY=angY(1:8000); subplot(2,1,1);plot(f,magY) xlabel(f);ylabel(|Y|);title(频域幅度谱) subplot(2,1,2);plot(f,angY) xlabel(f);ylabel(pi);title(频域相位谱)由绘图命令可得频域幅度谱和相位谱，如图4.2-3所示图4.2-3 频域幅度谱和相位谱4.3 设计椭圆滤波器语音信号解决时采用滤波器性能指标：fb=1100;fc=1200;As=20;Ap=1; fb=1100;fc=1200;As=20;Ap=1; wb=2*pi*fb/Fs;wc=2*pi*fc/Fs; wb=2*pi*fb/fs;wc=2*pi*fc/fs; T=1;OmegaP=(2/T)*tan(wb/2);OmegaS=(2/T)*tan(wc/2); c,d=afd_elip(OmegaP,OmegaS,Ap,As);* 椭圆滤波器阶次 = 4 b,a=bilinear(c,d,T); db,mag,pha,grd,w=freqz_m(b,a);由以上程序可得，双线性变换法设计椭圆低通滤波器幅度，幅度响应，群延时和相位响应等数据图，如图4.3所示图4.3 运用双线性变换设计椭圆低通滤波器（w（单位：）（a） 幅度（dB）；（b）幅度响应；（c）群延时；（d）相位响应4.4 信号滤波解决 signal=filter(b,a,y); subplot(2,1,1);plot(y) subplot(2,1,2);plot(signal)由以上程序可得，滤波先后时域波形对比图，如图4.4-1所示图4.4-1 滤波先后时域波形对比 X=fft(signal); magX=abs(X);angX=angle(X); subplot(2,2,1);plot(magX) subplot(2,2,2);plot(angX) subplot(2,2,3);plot(magY) subplot(2,2,4);plot(angY) sound (y,fs,bits)由以上程序可得，滤波先后频域频谱图，如图4.4-2所示所示图4.4-2 滤波先后频谱波形对比回放语音信号：将滤波后语音回放：sound(y,fs,bits)播放时没有听到含尖锐单频啸叫声，阐明已达到设计目。心得体会整个设计过程，让我感触最深就是功能强大性与掌握编程各种函数和语句重要性。固然一方面要理解所要编程运营对象原理。在课程设计过程中，我深深感受到我所学东西太少了，需要学习东西太多了。在一周课程设计时间里，我每天都过很充实，查资料、读程序，重复揣摩，学习过程是艰辛，但是同步也是高兴。通过实实在在课程设计，发现自己在课堂上所学知识对于解决实际问题来说，是远远不够，要想掌握技术，还需要加倍努力。但是，在这短短一周时间里，还是让我对数字信号解决有了更深理解，也学到了新知识让我掌握了visio、MATLAB等软件简朴使用，明白了软件仿真对设计重要性。再次，这次课程设计让我充分结识到团队合伙重要性，只有分工协作才干保证整个项目有条不紊。在整个课程设计过程中我再次结识到坚持、耐心、细心等品质重要，这对此后学习和工作是有很大协助。总之，这次课程设计让我获益良多。参照文献【1】matlab信号解决详解 陈亚勇等编著 人民邮电出版社，【2】电子滤波器设计 宁彦卿等译 科学出版社. 【3】 数字信号解决教程（第四版） 程佩青编著 清华大学出版社 【4】 数字椭圆滤波器Matlab设计与实现 王靖 李永全 当代电子技术第6 期总第245附录附录一：椭圆低通滤波器设计程序：%原始混合信号产生及对其采样Fs=100;t=(1:100)/Fs;s1=sin(2*pi*t*5);s2=sin(2*pi*t*15);s3=sin(2*pi*t*30);s=s1+s2+s3;figure(1);subplot(2，1，1)plot(t，s)xlabel(时间/t)ylabel(幅值)title(原始输入模仿信号)subplot(2，1，2)stem(t，s)xlabel(时间/t)ylabel(幅值)title(采样后输入信号)%椭圆低通滤波器设计%求取最小阶数和截止频率Wp=5;Ws=10;Fs=100;rp=0.1;rs=40;wp=2*Wp/Fs;ws=2*Ws/Fs;n，wn=ellipord(wp，ws，rp，rs)b，a=ellip(4，0.1，40，5*2/Fs);H，w=freqz(b，a，512);figure(2);plot(w*Fs/(2*pi)，abs(H);xlabel(频率/Hz);ylabel(幅度);title(椭圆低通滤波器幅度响应)grid;%对滤波后信号进行分析和变换sf=filter(b，a，s);figure(3);plot(t，sf);xlabel(时间/t);ylabel(幅值);title(滤波后信号)axis(0 1 -1 1);S=fft(s，512);SF=fft(sf，512);w=(0:255)/256*(Fs/2);figure(4);subplot(2，1，1)plot(w，abs(S(1:256);xlabel(频率/Hz);ylabel(幅度);title(滤波前傅立叶变换图)grid;subplot(2，1，2)plot(w，abs(SF(1:256);xlabel(频率/Hz);ylabel(幅度);title(滤波后傅立叶变换图)grid;附录二：第四章源程序清单 y,fs,bits = wavread(hy.wav); sound (y,fs,bits) plot(y);title(时域波形) t=(1:16000)/8000; plot(t,y);xlabel(t);ylabel(y);title(时域波形) Y=fft(y); magY=abs(Y);angY=angle(Y); w=(1:16000)/16000*2*pi; plot(w/pi,magY) f=w/(2*pi)*Fs; f=f(1:8000);magY=magY(1:8000);angY=angY(1:8000); subplot(2,1,1);plot(f,magY) xlabel(f);ylabel(|Y|);title(频域幅度谱) subplot(2,1,2);plot(f,angY) xlabel(f);ylabel(pi);title(频域相位谱) fb=1100;fc=1200;As=20;Ap=1; wb=2*pi*fb/Fs;wc=2*pi*fc/Fs; wb=2*pi*fb/fs;wc=2*pi*fc/fs; T=1;OmegaP=(2/T)*tan(wb/2);OmegaS=(2/T)*tan(wc/2); c,d=afd_elip(OmegaP,OmegaS,Ap,As);* 椭圆滤波器阶次 = 4 b,a=bilinear(c,d,T); db,mag,pha,grd,w=freqz_m(b,a); wc/pians = 0.3000 subplot(2,2,1);plot(w/pi,db);grid axis(0,1,-50,2); xlabel(w);ylabel(dB); subplot(2,2,2);plot(w/pi,mag);grid xlabel(w);ylabel(|Y|); subplot(2,2,3);plot(w/pi,grd);grid xlabel(w);ylabel(样本); subplot(2,2,4);plot(w/pi,pha);grid xlabel(w);ylabel(相位（单位：pi）); signal=filter(b,a,y); subplot(2,1,1);plot(y) subplot(2,1,2);plot(signal) X=fft(signal); magX=abs(X);angX=angle(X); subplot(2,2,1);plot(magX) subplot(2,2,2);plot(angX) subplot(2,2,3);plot(magY) subplot(2,2,4);plot(angY) sound (y,fs,bits)