数字信号处理IIR滤波器研究设计.doc

课 程 设 计 报 告课程名称数字信号处理课题名称IIR滤波器设计专 业 通信工程班 级 通信工程学 号姓 名指导教师胡瑛 张细政2013年 9月8日湖南工程学院课 程 设 计 任 务 书课程名称数字信号处理课 题IIR滤波器设计专业班级 通信工程学生姓名 学 号 指导老师 胡瑛 张细政审 批 任务书下达日期 2013 年 9月 1日任务完成日期 2013年 9月 8日数字信号处理课程设计任务书一 、设计目的综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用 MATLAB 作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。二、设计要求1、 MATLAB 的使用，掌握 MATLAB 的程序设计方法。2、 Windows 环境下语音信号采集的方法。3、数字信号处理的基本概念、基本理论和基本方法。4、TLAB 设计 FIR 和 nR 数字滤波器的方法。5、用 MATLAB 对信号进行分析和处理6、计报告4000以上，含程序设计说明，用户使用说明，源程序清单及程序框图。7、机演示。8、有详细的文档。文档中包括设计思路、设计仿真程序、仿真结果及相应的分析与结论。三、进度安排第一周星期一：课题讲解，查阅资料星期二: 总体设计，详细设计星期三：编程，上机调试、修改程序星期四：上机调试、完善程序星期五：答辩星期六-星期天：撰写课程设计报告附：课程设计报告装订顺序：封面、任务书、目录、正文、评分、附件（A4大小的图纸及程序清单）。正文的格式:一级标题用3号黑体,二级标题用四号宋体加粗,正文用小四号宋体;行距为22。正文的内容:一、课题的主要功能；二、课题的功能模块的划分；三、主要功能的实现；四、程序调试；五、总结；六、附件（所有程序的原代码，要求对程序写出必要的注释）；七、评分表。目录1、课题的主要功能11.1问题描述11.2功能要求12.、课题的功能模块的划分22.1语音信号采集（Butterworth型低通滤波器为例）22.2 语音信号分析22.3 含噪语音信号合成22.4 数字滤波器设计及滤波22.5回放语音信号33、滤波器实现的原理34、程序调试54.1 Butterworth型低通滤波器54.2Butterworth型高通滤波器64.3切比雪夫I型低通滤波器74.4切比雪夫I型带通滤波器84.5切比雪夫型带通滤波器94.6切比雪夫型低通滤波器104.7椭圆型的低通114.8 椭圆型的高通125、总结136、附件147、评分表241、课题的主要功能1.1问题描述采用IIR滤波器的完全设计函数来设计各型IIR滤波器（Butterworth型与切比雪夫I型、切比雪夫型和椭圆型的低通、高通、带通、带阻中的至少2种类型）,绘出滤波器的频域响应；并用这些数字滤波器对含噪语音信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。1.2功能要求1、语音信号采集（每个同学必须录制本人的语音信号）录制一段课程设计学生的语音信号并保存为文件，要求长度不小于10秒，并对录制的信号进行采样；录制时可以使用Windows自带的录音机，或者使用其它专业的录音软件，录制时需要配备录音硬件（如麦克风），为便于比较，需要在安静、干扰小的环境下录音。2、语音信号分析使用MATLAB绘出采样后的语音信号的时域波形和频谱图。根据频谱图求出其带宽，并说明语音信号的采样频率不能低于多少赫兹。3、含噪语音信号合成在MATLAB软件平台下，给原始的语音信号叠加上噪声，噪声类型分为如下几种：（1）白噪声；（2）单频噪色（正弦干扰）；（3）多频噪声（多正弦干扰）；（4）其它干扰，可设置为低频、高频、带限噪声，或Chirp干扰、冲激干扰。绘出叠加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析并体会含噪语音信号频谱和时域波形的改变。4、数字滤波器设计及滤波利用不同的原型低通滤波器（Butterworth型与切比雪夫I型、切比雪夫型和椭圆型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪语音信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。5、回放语音信号对滤波后的语音信号进行回放，感觉滤波前后语音信号的变化。2.、课题的功能模块的划分2.1语音信号采集（Butterworth型低通滤波器为例）x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);2.2语音信号分析原始信号x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)2.3含噪语音信号合成加噪后的语音信号N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)2.4数字滤波器设计及滤波加噪后的语音信号进行滤波处理Fs=fs; %采样频率wp=4000*2/Fs;ws=4300*2/Fs; %根据采样频率将滤波器边界频率进行转换Rp=1;Rs=20; %通带波纹和阻带衰减Nn=128; %显示滤波器频率特性的数据长度N,Wn=buttord(wp,ws,Rp,Rs); %求得数字滤波器的最小阶数和截止频率（归一化频率）b,a=butter(N,Wn); %设计Butterworth高通数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filter(b,a,x3); %对输入信号进行滤波2.5回放语音信号subplot(2,1,1)plot(y);title(滤波语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(滤波语音信号频谱)%sound(y,fs,bits);3、滤波器实现的原理IIR数字滤波器设计原理：IIR数字滤波器是通过因果稳定的Ha（s）映射成因果稳定的H（z），即s平面的左半平面必须映射到z平面单位圆的内部。H（z）的频率响应能模仿Ha（s）的频率响应，即s平面的虚轴必须映射到z平面的单位圆上。变换前后的滤波器在时域或频域的主要特征（频率响应或单位冲激响应等）应尽可能相同或接近。将传输函数Ha（s）从s平面转换到z平面的方法有多种，主要有冲激不变法和双线性变换法。如果不考虑混叠现象，利用冲击不变法实现的数字滤波器会很好地重现原模拟滤波器的频率特性。而且数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域特性逼近好。在本设计中，用冲激函数作为系统激励信号，用各种数字滤波器作为测试系统。冲激函数具有无限宽广的频谱，用冲激函数做激励信号相当于对测试系统输入所有频率的信号，系统必然有对应的输出。用Transfer函数计算出系统输出与输入的傅立叶变换之比，从而得到系统的频率响应函数。IIR数字滤波器的系统函数为的有理分式：设计IIR滤波器的系统函数，就是要确定的阶数N及分子分母多项式的系数和，使其满足指定的频率特性。由于模拟滤波器的设计有许多简单而严谨的设计公式和大量的图表可以利用，因此IIR滤波器设计的方法之一是：先设计一个合适的模拟滤波器，然后将模拟滤波器通过适当的变换转换成满足给定指标的数字滤波器。1、Butterworth模拟低通滤波器幅度平方函数：其中，N为滤波器的阶数，为通带截止频率。2Chebyshev模拟低通滤波器IIR数字滤波器经典设计法的一般步骤为：根据给定的性能指标和方法不同，首先对设计性能指标中的频率指标，如数字边界频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。估计模拟滤波器最小阶数和截止频率，利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。设计模拟低通滤波器原型。利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等），利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bilinear或impinvar。后面的几个步骤在前面的章节中已经论述，这里主要介绍第一步，关于设计性能指标的转换。设计IIR滤波器时，给出的性能指标通常分数字指标和模拟指标两种。数字性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。数字频率和的取值范围为0，单位弧度。而MATLAB工具函数常采用归一化频率，和的取值范围为01，对应于0，此时需进行转换。模拟性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。模拟频率和单位为弧度/秒（rad/s）。MATLAB信号处理工具箱中，设计性能指标的转换应根据不同设计方法进行不同处理。4、程序调试4.1 Butterworth型低通滤波器 图4-1-1 图4-1-2图4-1-3 图4-1-44.2 Butterworth型高通滤波器图4-2-1 图4-2-2图4-2-3 图4-2-44.3 切比雪夫I型低通滤波器 图4-3-1 图4-3-2图4-3-3 图4-3-44.4 切比雪夫I型带通滤波器图4-4-1 图4-4-2图4-4-3 图4-4-44.5 切比雪夫型带通滤波器 图4-5-1 图4-5-2 图4-5-3 图4-5-44.6 切比雪夫型低通滤波器 图4-6-1 图4-6-2 图4-6-3 图4-6-44.7 椭圆型的低通图4-7-1 图4-7-2图4-7-3 图4-7-44.8 椭圆型的高通图4-8-1 图4-8-2图4-8-3 图4-8-45、总结通过这次课程设计，我对matlab软件更加了解，更加熟练，用得更加顺手，对数字信号中的滤波器也有了新的认识，本次我是采用IIR滤波器的完全设计函数来设计各型IIR滤波器（Butterworth型与切比雪夫I型、切比雪夫型和椭圆型的低通、高通、带通、带阻中的至少2种类型）,并要求绘出滤波器的频域响应；并用这些数字滤波器对含噪语音信号分别进行滤波处理，它的设计方法：（1）根据给定的性能指标和方法不同，首先对设计性能指标中的频率指标，如数字边界频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。（2）估计模拟滤波器最小阶数和截止频率，利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。(3)设计模拟低通滤波器原型。利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。(4)由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等），利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。(5)将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bilinear或impinvar。数字滤波器可以通过编程实现各种不同系统，满足不同的需要，也可以随时改动系数，调整滤波器参数，选择最佳方案。使用LabVIEW软件平台开发电气参数测量仪等虚拟仪器，实现了更高的效率，节省了更多的硬件开销，方便了系统的维护和减轻了仪器更新的负担。使用虚拟仪器逐步代替传统仪器已经成为测试领域发展的趋势。但是在实际应用中，仍要根据具体情况进行程序的优化和软硬件的结合，使虚拟仪器发挥更高的性能。6、附件Butterworth型低通滤波器x1,fs,bits=wavread(d:/baojing2);%sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=4000*2/Fs;ws=4300*2/Fs; %根据采样频率将滤波器边界频率进行转换Rp=1;Rs=20; %通带波纹和阻带衰减Nn=128; %显示滤波器频率特性的数据长度N,Wn=buttord(wp,ws,Rp,Rs); %求得数字滤波器的最小阶数和截止频率（归一化频率）b,a=butter(N,Wn); %设计Butterworth高通数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filter(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(滤波语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(滤波语音信号频谱)%sound(y,fs,bits);Butterworth型高通滤波器x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*40*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=1000; %采样频率wp=50*2/Fs;ws=40*2/Fs; %根据采样频率将滤波器边界频率进行转换Rp=1;Rs=20; %通带波纹和阻带衰减Nn=128; %显示滤波器频率特性的数据长度N,Wn=buttord(wp,ws,Rp,Rs); %求得数字滤波器的最小阶数和截止频率（归一化频率）b,a=butter(N,Wn,high); %设计Butterworth高通数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filter(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(原始语音信号频谱)%sound(y,fs,bits);切比雪夫I型带通滤波器x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=50 4000*2/Fs; %通带边界频率（归一化频率）（6-20式）ws=30 5000*2/Fs; %阻带边界频率（归一化频率）（6-20式）Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数N,Wn=cheb1ord(wp,ws,Rp,Rs);%求得数字滤波器的最小阶数和归一化截止频率b,a=cheby1(N,Rp,Wn); %按最小阶数、通带波纹和截止频率设计数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %求得滤波器的频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filtfilt(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(滤波语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(滤波语音信号频谱)%sound(y,fs,bits);切比雪夫I型低通滤波器x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=4000*2/Fs; %通带边界频率（归一化频率）（6-20式）ws=4500*2/Fs; %阻带边界频率（归一化频率）（6-20式）Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数N,Wn=cheb1ord(wp,ws,Rp,Rs);%求得数字滤波器的最小阶数和归一化截止频率b,a=cheby1(N,Rp,Wn); %按最小阶数、通带波纹和截止频率设计数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %求得滤波器的频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filtfilt(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(原始语音信号频谱)%sound(y,fs,bits);切比雪夫型低通滤波器x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=4000*2/Fs; ws=4500*2/Fs; Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数N,Wn=cheb1ord(wp,ws,Rp,Rs);%求得数字滤波器的最小阶数和归一化截止频率b,a=cheby1(N,Rp,Wn); %按最小阶数、通带波纹和截止频率设计数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %求得滤波器的频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filtfilt(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(滤波语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(滤波语音信号频谱)%sound(y,fs,bits);切比雪夫型带通滤波器x1,fs,bits=wavread(d:/baojing2);sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=50 4000*2/Fs; %通带边界频率（归一化频率）（6-20式）ws=40 5000*2/Fs; %阻带边界频率（归一化频率）（6-20式）Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数N,Wn=cheb2ord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率b,a=cheby2(N,Rs,Wn); %按最小阶数截止频率和阻带衰减设计数字滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %按传递函数系数、数据点数和采样频率求频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filtfilt(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(原始语音信号频谱)%sound(y,fs,bits);椭圆型的高通x1,fs,bits=wavread(d:/baojing2);%sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*40*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);axis(0,35000,0,400)title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=50*2/Fs;ws=40*2/Fs; %根据采样频率将滤波器边界频率进行转换Rp=1;Rs=20; %通带波纹和阻带衰减Nn=128; %显示滤波器频率特性的数据长度N,Wn=ellipord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率b,a=ellip(N,Rp,Rs,Wn,high); %按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filtfilt(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(滤波语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(滤波语音信号频谱)%sound(y,fs,bits);椭圆型的低通x1,fs,bits=wavread(d:/baojing2);%sound(x1,fs,bits);y1=fft(x1,32768);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y1);title(原始语音信号频谱)N=length(x1)-1;t=0:1/fs:N/fs;d=0.8*cos(2*pi*5000*t);x3=x1+d;%sound(x3,fs,bits);figure(2)subplot(2,1,1)plot(x3)title(加单频余弦信号的语音信号时域波形)y3=fft(x3,32768);subplot(2,1,2)plot(abs(y3);title(加单频余弦信号的语音信号频谱)Fs=fs; %采样频率wp=4000*2/Fs;ws=5000*2/Fs; %根据采样频率将滤波器边界频率进行转换Rp=1;Rs=20; %通带波纹和阻带衰减Nn=128; %显示滤波器频率特性的数据长度N,Wn=ellipord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率b,a=ellip(N,Rp,Rs,Wn); %按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器figure(3)H,f=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(4)y=filtfilt(b,a,x3); %对输入信号进行滤波y4=fft(y,32768);subplot(2,1,1)plot(y);title(原始语音信号时域波形);subplot(2,1,2)plot(abs(y4);title(原始语音信号频谱)%sound(y,fs,bits);Butterworth型低通滤波器和Butterworth型