MATLAB课程设计（论文）FIR滤波器的分析及设计.doc

1 引言1.1 MATLAB的概况MATLAB是矩阵实验室（Matrix Laboratory）之意。除具备优越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建筑仿真和实时控制等功能。MATLAB的基本数据单位是矩阵，它的指令表达式与数学，工程中常用的形式十分相似，故用MATLAB来解算问题要比C、FORTRAN等语言相同的事情简捷得多。MATLAB包括拥有数百个内部函数的主包和三十几种工具包（Toolbox）。工具包又可以分为功能性工具包和科学工具包。开放性使MATLAB广受用户欢迎。除内部函数外，所有MATLAB主包文件和各种工具包都是可读可修改的文件，用户通过对源程序的修改或加入自己编写程序构造新的专用工具包。1.2 MATLAB的语言特点MATLAB最突出的特点就是简洁。它用更直观的、符合人们思维习惯的代码，代替了C和FORTRAN语言的冗长代码。它给用户带来的是最直观，最简洁的程序开发环境。以下简介MATLAB的主要特点：（1） 语言简洁紧凑，使用方便灵活，库函数及其丰富。程序书写形式自由，利用丰富的库函数避开繁杂的子程序编写任务，压缩了一切不必要的编程工作。（2） 运算符丰富，由于MATLAB使用C语言编写的，它提供了和C语言几乎一样多的运算符，灵活使用MATLAB的运算符将使程序变得极为简短。（3） MATLAB既具有结构化的控制语句（如for循环，while循环，break语句和if语句），又有面向对象编程的特性。（4） 程序限制不严格，程序设计自由度大。例如，在MATLAB里，用户无需对矩阵预定义就可使用。（5） 程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。（6） MATLAB的图形功能强大，在FORTRAN和C语言里，绘图都很不容易，但在MATLAB里，数据的可视化非常简单。它还具有较强的编辑图形界面的能力。2 FIR滤波器的分析及设计2.1 用窗函数法设计FIR滤波器的步骤如下：(1)根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数）。窗函数类型可根据最小阻带衰减AS独立选择，因为窗口长度N对最小阻带衰减AS没有影响。在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N。设待求滤波器的过渡带宽为，它与窗口长度N近似成反比。窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正。原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N。在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数。(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应。如果给出待求滤波器的频率响应为，则理想的单位脉冲响应可以用傅里叶反变换式得到，在一般情况下，是不能用封闭公式表示的，需要采用数值方法表示。从到采样N点，采用离散傅里叶反变换(IDFT)即可求出。(3)计算滤波器的单位脉冲响应h(n)。它是理想单位脉冲响应和窗函数的乘积，即，在MATLAB中用点乘命令表示为。(4)验算技术指标是否满足要求。为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。2.2 FIR滤波器的设计2.2.1 FIR低通滤波器详细设计过程 (1) 确定FIR模拟低通滤波器的技术指标： (2) 归一化，求得数字边缘频率： (3) 计算中心频率、过度带宽，设定旁带宽度： (4) 计算beta和N： (5) 利用fir1函数确定凯泽窗低通滤波器的实际单位脉冲响应： (6) 利用freqz函数求得幅频特性： (7) 作出的幅频特性曲线： (8) 利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例采用的是一段音乐，采样频率为22050Hz： (9) 利用sound函数对原始语音信号进行回放： (10) 利用filter函数对样本信号进行滤波处理： (11) 利用sound函数对经过滤波器后的语音信号进行回放： (12) 利用plot函数分别作出信号在滤波前后的时域波形： (13) 利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换： (14) 利用plot函数分别作出信号在滤波前后的频域波形： 2.2.2程序流程图，如图1所示：开始确定模拟低通滤波器的设计指标技术指标从模拟转换到数字计算中心频率，用窗函数法设计FIR滤波器的步骤如下：(1)根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1）。窗函数类型可根据最小阻带衰减AS独立选择，因为窗口长度N对最小阻带衰减AS没有影响。在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N。设待求滤波器的过渡带宽为，它与窗口长度N近似成反比。窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正。原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N。在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数。(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应。如果给出待求滤波器的频率响应为，则理想的单位脉冲响应可以用傅里叶反变换式得到，在一般情况下，是不能用封闭公式表示的，需要采用数值方法表示。从=0到=2采样N点，采用离散傅里叶反变换(IDFT)即可求出。(3)计算滤波器的单位脉冲响应h(n)。它是理想单位脉冲响应和窗函数的乘积，即，在MATLAB中用点乘命令表示为。(4)验算技术指标是否满足要求。为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。过渡带宽，alfa用fir1函数确定滤波器的实际脉冲响应h(n)用freqz确定幅频响应做出凯泽窗的幅频响应曲线结束图1 FIR低通滤波器的程序流图2.2.3凯泽窗低通滤波器的幅频特性，如图2所示：图2 FIR低通滤波器的幅频特性 2.2.4输入信号、输出信号的时域波形如图3、图4所示： 图3 输入信号的时域波形图4 输出信号的时域波形 2.2.5语音信号滤波前后的频谱图，如图5、图6所示：图5 语音信号滤波前的频谱图图6 语音信号滤波后的频谱图分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的低频率段被保留，高频率段被滤除。2.2.6 FIR带通滤波器详细设计过程 (1) 确定FIR模拟带通滤波器的技术指标： (2) 归一化，求得数字边缘频率： (3) 计算过度带宽，设定旁带宽度： (4) 计算beta和N： (5) 利用fir1函数确定凯泽窗带通滤波器的实际单位脉冲响应： (6) 利用freqz函数求得幅频特性： (7) 作出的幅频特性曲线： (8) 利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例采用的是一段音乐，采样频率为22050Hz： (9) 利用sound函数对原始语音信号进行回放： (10) 利用filter函数对样本信号进行滤波处理： (11) 利用sound函数对经过滤波器后的语音信号进行回放： (12) 利用plot函数分别作出信号在滤波前后的时域波形： (13) 利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换： (14) 利用plot函数分别作出信号在滤波前后的频域波形： 2.2.7程序流程图，如图7所示：图7 FIR带通滤波器设计的程序流图开始确定模拟带通滤波器的设计指标技术指标从模拟转换到数字计算过渡带宽，alfa，用窗函数法设计FIR滤波器的步骤如下：(1)根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1）。窗函数类型可根据最小阻带衰减AS独立选择，因为窗口长度N对最小阻带衰减AS没有影响。在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N。设待求滤波器的过渡带宽为，它与窗口长度N近似成反比。窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正。原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N。在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数。(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应。如果给出待求滤波器的频率响应为，则理想的单位脉冲响应可以用傅里叶反变换式得到，在一般情况下，是不能用封闭公式表示的，需要采用数值方法表示。从=0到=2采样N点，采用离散傅里叶反变换(IDFT)即可求出。(3)计算滤波器的单位脉冲响应h(n)。它是理想单位脉冲响应和窗函数的乘积，即，在MATLAB中用点乘命令表示为。(4)验算技术指标是否满足要求。为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。过渡带宽，alfa用fir1函数确定滤波器的实际脉冲响应h(n)用freqz确定幅频响应做出凯泽窗的幅频响应曲线结束2.2.8凯泽窗带通滤波器的幅频特性，如图8所示：图8 FIR带通滤波器的幅频特性 2.2.9输入信号、输出信号的时域波形如图9、图10所示：图9 输入信号的时域波形图10 输出信号的时域波形2.2.10语音信号滤波前后的频谱图，如图11、图12所示：图11 语音信号滤波前的频谱图图12 语音信号滤波后的频谱图分析滤波前后信号的变化，波形图变窄。频谱变化表明:信号的低频率段和高频率段被滤除，中间频率段被保留。2.2.11 FIR高通滤波器详细设计过程 (1) 确定FIR模拟高通滤波器的技术指标： (2) 归一化，求得数字边缘频率： (3) 计算过度带宽，设定旁带宽度： (4) 计算beta和N： (5) 利用fir1函数确定凯泽窗高通滤波器的实际单位脉冲响应： (6) 利用freqz函数求得幅频特性： (7) 作出的幅频特性曲线： (8) 利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例采用的是一段音乐，采样频率为22050Hz： (9) 利用sound函数对原始语音信号进行回放： (10) 利用filter函数对样本信号进行滤波处理： (11) 利用sound函数对经过滤波器后的语音信号进行回放： (12) 利用plot函数分别作出信号在滤波前后的时域波形： (13) 利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换： (14) 利用plot函数分别作出信号在滤波前后的频域波形： 2.2.12程序流程图，如图13所示：开始确定模拟高通滤波器的设计指标技术指标从模拟转换到数字计算中心频率，过渡带宽，alfa用fir1函数确定滤波器的实际脉冲响应h(n)用freqz确定幅频响应做出凯泽窗的幅频响应曲线结束图13 FIR高通滤波器设计的程序流图2.2.13凯泽窗高通滤波器的幅频特性，如图14所示：图14 FIR高通滤波器的幅频特性2.2.14输入信号、输出信号的时域波形如图15、图16所示：图15 输入信号的时域波形图16 输出信号的时域波形2.2.15语音信号滤波前后的频谱图，如图17、图18所示：图17 语音信号滤波前的频谱图图18 语音信号滤波后的频谱图分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的高频率段被保留，低频率段被滤除。3 IIR滤波器的分析及设计3.1 用双线性变换法设计IIR滤波器的步骤如下： 利用模拟滤波器来设计IIR数字滤波器，就是从已知的模拟滤波器系统函数设计IIR数字滤波器系统函数，也就是从S平面到Z平面的变换，这个变换通常是复变函数的映射变换，这个变换需要满足如下两个基本要求：（1）的频响应该模拟的频响，即将S平面的虚轴映射到Z平面的单位圆上。（2）因果稳定的通过映射后的仍应该因果稳定，即S平面的左半平面即应该被映射到Z平面单位圆以内，即|Z|<1。由变换成有两种方法：脉冲响应不变法和双线性变换法。本次课程设计用双线性变换法来设计IIR滤波器。双线性变换法可以实现由S平面一对一的映射到Z平面，它通过一个过渡实现，过渡的面设为，从S平面到平面的映射是将频率非线性压缩，将S平面整个频率轴上的频率范围压缩到平面的之间，从平面到Z平面的映射是采用的转换关系变换的，它是一对一的进行转换，从而避免了混叠现象。S平面的虚轴上，平面的虚轴上，为使S平面中的由经过0变化到时，平面的从经过零变化到，采用的频率变换关系式。从而实现了S平面上整个虚轴完全压缩到平面上虚轴的之间推广到S平面到平面的转换为。而平面仍然通过关系式转换到Z平面，从S平面平面Z平面，得到S与Z的关系为或从而实现了S平面到Z平面的转换，即双线性变换法。3.2 IIR滤波器的设计3.2.1 IIR低通滤波器的详细设计过程 (1) 确定IIR模拟低通滤波器的技术指标： (2) 归一化，求得数字边缘频率： (3) 技术指标从数字到模拟转换： (4) 利用buttord函数求得低通滤波器的阶数N和3dB截止频率wc： (5) 利用buttap函数求得归一化模拟低通滤波器的零点，极点和增益： (6) 利用zp2tp函数求得模拟低通滤波器的系统函数： (7) 利用lp2lp函数，令wc=wap,反归一化： (8)利用bilinear函数实现到的转化： (9)利用freqz函数求得的幅频特性： (10)作出的幅频特性曲线： (11) 利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例采用的是一段音乐，采样频率为22050Hz： (12) 利用sound函数对原始语音信号进行回放： (13) 利用filter函数对样本信号进行滤波处理： (14) 利用sound函数对经过滤波器后的语音信号进行回放： (15) 利用plot函数分别作出信号在滤波前后的时域波形： (16) 利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换，采用一次基2按时间抽取： (14) 利用plot函数分别作出信号在滤波前后的频域波形： 3.2.2 程序流程图，如图19所示： 结束开始确定模拟低通滤波器的技术指标归一化，求数字边缘频率技术指标从数字到模拟确定滤波器的阶数N和wc计算归一化频率的零点、极点和增益模拟低通滤波器的系统函数利用lp2lp函数，令wc=wap,反归一化双线性变换法，从s域到z域H（z）的幅频特性做出H（z）的幅频特性曲线图19 巴特沃兹低通滤波器设计的程序流图3.2.3 巴特沃兹低通滤波器的幅频特性，如图20所示：图20 巴特沃兹低通滤波器的幅频特性3.2.4信号滤波前后的时域波形，如图21、22所示：图21 滤波前的时域波形图22 滤波后的时域波形3.2.5 原语音信号滤波前后的频谱图，如图23、24所示：图23 原语音信号滤波前的频谱 图24 原语音信号经巴特沃兹低通滤波器滤波前的频谱分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的低频率段被保留，高频率段被滤除。3.2.6 IIR带通滤波器的详细设计过程 (1) 确定IIR模拟带通滤波器的技术指标： (2) 归一化，求得数字边缘频率： (3) 技术指标从数字到模拟转换： (4) 利用buttord函数求得带通滤波器的阶数N和3dB截止频率wc： (5) 利用buttap函数求得归一化模拟带通滤波器的零点，极点和增益： (6) 利用zp2tp函数求得模拟带通滤波器的系统函数： (7) 利用lp2bp函数，令wc=wap,反归一化： (8)利用bilinear函数实现到的转化： (9)利用freqz函数求得的幅频特性： (10)作出的幅频特性曲线： (11) 利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例采用的是一段音乐，采样频率为22050Hz： (12) 利用sound函数对原始语音信号进行回放： (13) 利用filter函数对样本信号进行滤波处理： (14) 利用sound函数对经过滤波器后的语音信号进行回放： (15) 利用plot函数分别作出信号在滤波前后的时域波形： (16) 利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换，采用一次基2按时间抽取： (14) 利用plot函数分别作出信号在滤波前后的频域波形： 3.2.7 程序流程图，如图25所示：图25 IIR带通滤波器设计的程序流图开始确定模拟带通滤波器的技术指标归一化，求数字边缘频率技术指标从数字到模拟确定滤波器的阶数N和wc计算归一化频率的零点、极点和增益模拟带通滤波器的系统函数利用lp2bp函数，令wc=wap,反归一化双线性变换法，从s域到z域H（z）的幅频特性做出H（z）的幅频特性曲线结束3.2.8 巴特沃兹带通滤波器的幅频特性，如图26所示：图26 巴特沃兹带通滤波器的幅频特性3.2.9信号滤波前后的时域波形，如图27、28所示：图27 信号滤波前的时域波形图28 信号滤波后的时域波形3.2.10 原语音信号滤波前后的频谱图，如图29、30所示：图29 信号滤波前的频谱图图30 信号滤波后的频谱图分析滤波前后信号的变化，波形图变窄。频谱变化表明:信号的低频率段和高频率段被滤除，中间频率段被保留。3.2.11 IIR高通滤波器的详细设计过程 (1) 确定IIR模拟高通滤波器的技术指标： (2) 归一化，求得数字边缘频率： (3) 技术指标从数字到模拟转换： (4) 利用buttord函数求得高通滤波器的阶数N和3dB截止频率wc： (5) 利用buttap函数求得归一化模拟高通滤波器的零点，极点和增益： (6) 利用zp2tp函数求得模拟高通滤波器的系统函数： (7) 利用lp2lp函数，令wc=wap,反归一化： (8)利用bilinear函数实现到的转化： (9)利用freqz函数求得的幅频特性： (10)作出的幅频特性曲线： (11) 利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例采用的是一段音乐，采样频率为22050Hz： (12) 利用sound函数对原始语音信号进行回放： (13) 利用filter函数对样本信号进行滤波处理： (14) 利用sound函数对经过滤波器后的语音信号进行回放： (15) 利用plot函数分别作出信号在滤波前后的时域波形： (16) 利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换，采用一次基2按时间抽取： (14) 利用plot函数分别作出信号在滤波前后的频域波形： 3.2.12 程序流程图，如图31所示：开始确定模拟高通滤波器的技术指标归一化，求数字边缘频率技术指标从数字到模拟确定滤波器的阶数N和wc计算归一化频率的零点、极点和增益模拟高通滤波器的系统函数利用lp2hp函数，令wc=wap,反归一化双线性变换法，从s域到z域H（z）的幅频特性做出H（z）的幅频特性曲线结束图31 IIR高通滤波器设计的程序流图3.2.13 巴特沃兹高通滤波器的幅频特性，如图32所示：图32 巴特沃兹高通滤波器的幅频特性3.2.14信号滤波前后的时域波形，如图33、34所示：图33 信号滤波前的时域波形图34 信号滤波后的时域波形3.2.15 原语音信号滤波前后的频谱图，如图35、36所示：图35 原语音信号滤波前的频谱图图36 原语音信号滤波后的频谱图分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的高频率段被保留，低频率段被滤除。4 图形用户界面GUI图形用户界面（Graphical User Interfaces,GUI）是一种提供人机交互的工具和方法。GUI是包含图形对象，如窗口、图标、菜单和文本等图文并茂的用户界面。4.1 启动GUIDE，界面如图37所示：图37 GUIDE的启动界面单击OK，新建一个GUI,如图38所示：图38 具有空白页面的GUIDE4.2 建立GUIDE，界面如图39所示：图39 建立好的GUIDE运行的GUI界面,如图40所示：图40 运行的GUIDE5 总结转眼间，两个星期过去了，课设也将结束了。想到刚刚拿到课设题目时，一头雾水，不知道该从哪里下手，到现在看着做出的成果，心中的喜悦之情油然而生。通过这两周的滤波器课程设计，更加深入了解了matlab的功能，以及各种滤波器的性能。MATLAB最突出的特点就是简洁。它用更直观的、符合人们思维习惯的代码，代替了C和FORTRAN语言的冗长代码。它给用户带来的是最直观，最简洁的程序开发环境。通过matlab软件将滤波器的波形更加形象直观的表示出来，使之更加容易理解，更加容易开发研究。但是仅仅两周的时间是远远不够熟练掌握MATLAB软件的使用的，因此，在今后的时间里，我会更加努力的学习matlab，以便可以做出一些成就，制作一些作品，使所学知识能够熟练掌握并将其运用到实际中来。最后，我要感谢老师和同学的帮助，我现在的成果和你们的帮助是密不可分的，我会更加努力的学习，不辜负大家的期望，谢谢！参考文献1范寿康.DSP技术与DSP芯片. 北京:电子工业出版社，20082万永革.数字信号处理的MATLAB实现. 北京:科学出版社， 20073程佩青.数字信号处理教程. 北京：清华大学出版社出版，20014赵健,李毅.数字信号处理学习辅导. 北京:清华大学出版社出版, 20065丁玉美.数字信号处理学习指导与题解. 北京:电子工业出版社出版, 20076丛玉良.数字信号处理原理及其MATLAB实现. 北京:电子工业出版社,20057陈生谭,郭宝龙等.信号与系统. 西安:西安电子科技大学出版社,20018王家文,王皓等.MATLAB7.0编程基础. 北京:机械工业出版社,20059李哲英,骆丽等.DSP基础理论与应用技术. 北京:航空航天大学出版社,2002附录%FIR低通滤波器的设计fp=1000;fs=1200; As=120;Fs=22050; wp=2*pi*fp/Fs; ws=2*pi*fs/Fs;wo=ws-wp; belta=0.1102*(As-0.87); N=ceil(As-8)/2.285/wo); wc=(wp+ws)/2/pi;b=fir1(N,wc,kaiser(N+1,belta);Hw,w=freqz(b,1);figure(1)plot(w*Fs/2/pi,20*log10(abs(Hw); grid;title('FIR低通滤波器');ylabel('幅度/db');xlabel('频率/Hz');axis(0,8000,-200,10);x=wavread('D:0072.wav'); sound(x,22050); y=filter(b,1,x); sound(y,22050); figure(2) subplot(2,1,1)plot(x(:,1); title('输入信号时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,1,2)plot(y(:,1); title('输出信号时域波形');xlabel('f/Hz');ylabel('幅度');figure(3)X=fft(x(:,1),1024); Y=fft(y(:,1),1024); f=Fs*0:1023/1024;subplot(2,1,1)plot(f,abs(X); title('原始信号的FFT波形');ylabel('幅度');xlabel('频率/Hz');axis(0,8000,0,4);subplot(2,1,2)plot(f,abs(Y);title('原始语音信号经滤波器滤波后的FFT');ylabel('幅度');xlabel('频率/Hz');axis(0,8000,0,4);%带通滤波器的设计fb1=1200;fc1=1000; fb2=3000;fc2=3200;As=120;Fs=22050;wp1=2*pi*fc1/Fs; ws1=2*pi*fb1/Fs;wp2=2*pi*fc2/Fs;ws2=2*pi*fb2/Fs;wo1=wp1-ws1; wo2=wp2-ws2;wo=max(wo1,wo2);beta=0.1102*(As-0.87); N=ceil(As-8)/2.285/wo); wn1=(fb1+fc1)/Fs;wn2=(fb2+fc2)/Fs;wn=wn1,wn2;b=fir1(N,wn,'band',kaiser(N+1,beta);Hw,w=freqz(b,1,512,Fs);figure(4);plot(w,20*log10(abs(Hw); grid;title('FIR带通滤波器');ylabel('幅度/db');xlabel('频率/Hz');axis(0,8000,-200,10);x=wavread('D:0072.wav'); sound(x,22050); y=filter(b,1,x); sound(y,22050); figure(5);subplot(2,1,1);plot(x(:,1); title('输入信号时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,1,2);plot(y(:,1);title('输出信号时域波形');xlabel('f/Hz');ylabel('幅度');figure(6);X=fft(x(:,1),1024); Y=fft(y(:,1),1024); f=Fs*0:1023/1024;subplot(2,1,1);plot(f,abs(X); title('原始语音信号的FFT波形');ylabel('幅度');xlabel('频率/Hz');axis(0,8000,0,3);subplot(2,1,2);plot(f,abs(Y); title('原始语音信号经滤波器滤波后的FFT波形');ylabel('幅度');xlabel('频率/Hz');axis(0,8000,0,3);title('输入信号时域波形');%高通滤波器的设计fp=5000;fs=4800;As=120;Fs=22050;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;wo=wp-ws;beta=0.1102*(As-0.87);N=ceil(As-8)/2.285/wo);wc=(wp+ws)/2/pi;b=fir1(N,wc,'high',kaiser(N+2,beta);Hw,w=freqz(b,1);figure(7);plot(w*Fs/2/pi,20*log10(abs(Hw);grid;title('FIR高通滤波器的幅频特性');ylabel('幅度/dB');xlabel('频率/Hz');axis(0,10000,-200,10);x=wavread('D:0072.wav');sound(x,22050);y=filter(b,1,x);sound(y,22050);figure(8);subplot(2,1,1);plot(x(:,1);title('输入信号的时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,1,2);plot(y(:,1);title('输出信号的时域波形');xlabel('f/Hz');ylabel('幅度');figure(9);X=fft(x(:,1),1024);Y=fft(y(:,1),1024);f=Fs*0:1023/1024;subplot(2,1,1);plot(f,abs(X);title('原始语音信号的FFT波形');ylabel('幅度');xlabel('频率/Hz');axis(0,10000,0,3);subplot(2,1,2);plot(f,abs(Y);title('原始语音信号经滤波器滤波后的FFT波形');ylabel('幅度');xlabel('频率/Hz');axis(0,10000,0,3);%IIR低通滤波器fp=1000;fs=1200;Fs=22050;Ap=1;As=15;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Ts=