语音信号滤波去噪——使用双线性变换法设计切比雪夫I型滤波器

1、语音信号滤波去噪双线性变换设计切比雪夫i型滤波器第21页 共21 页语音信号滤波去噪使用双线性变换法设计切比雪夫i型滤波器学生姓名：xxx 指导老师：xxx摘 要 本次课程设计的主要是设计一个切比雪夫i型滤波器，并对语音信号进行滤波去噪。此次设计用双线性变换法，开发平台为matlab集成系统。用windows下的录音机采集一段语音信号，绘制波形并观察其频谱。给定相应技术指标，设计一个满足指标的切比雪夫i型iir滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。程序通过调试运行，初步实现了设计目标，达到设计要求。关键词 课程设计；滤波去噪；iir滤波器；切比雪夫i型滤波器

2、；matlab1 引 言通过麦克风采集一段8000hz，8k的单声道语音信号，绘制其波形并观察频谱，利用双线性变换法，以给定的设计指标，通带截止频率1700hz，阻带截止频率1750hz，通带波纹1db，阻带波纹15db，设计一个切比雪夫i型iir滤波器，对采集的语音信号进行滤波去噪处理。1.1课程设计目的本次课程设计的目的是让学生综合数字信号处理和matlab并实现一个较为完整的小型滤波器系统。综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用 matlab 作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。通过系统设计、软件仿真、程序安排

3、与仿真、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2课程设计的要求（1）学会 matlab 的使用，掌握 matlab 的程序设计方法；（2）滤波器指标必须符合工程实际，根据模拟滤波器的性能指标，确定数字滤波器指标；（3）采用双线性变换法，设计满足上述性能指标要求的chebyshevi型数字低通滤波器；（4）设计完后应检查其频率响应曲线是否满足指标；（5）处理结果和分析结论应该一致，而且应符合理论；（6）独立完成课程设计并按要求编写课程设计报告书；1.3设计平台matlab是一种既可交互使用又能解释执行的计算机变成语言，可以用

4、直观的数学表达式来描述问题，从而避开繁琐的底层编程，并且是解决工程技术问题的计算平台。它拥有丰富的函数资源和工具箱资源，mtlab最大的特点就是简洁开放的程序代码和直观实用的开发环境。具备库函数资源丰富；语言精练，代码灵活；面向对象，控制功能优良；程序设计优良，图形功能强大；源代码开放，形形色色工具箱等优点。matlab已经成为国际上最流行的科学与工程计算软件工具1。2 设计原理2.1 iir滤波器从离散时间来看，若系统的单位抽样（冲激）响应延伸到无穷长，称之为“无限长单位冲激响应系统”，简称iir系统。无限长单位冲激响应滤波器有以下及个特点：（1）系统的单位冲激响应h(n)是无限长；（2）系

5、统函数h(z)在有限z平面(0|z|1时, |x|, vn(x) 切比雪夫滤波器的振幅平方特性如图2.1所示通带内，振幅平方函数的变化范围为n为偶数 （2-4） n为奇数 （2-5） 切比雪夫滤波器的振幅平方特性如图1所示图2.1 切比雪夫滤波器的振幅平方特性图 2.3 双线性变换法双线性变换由下式定义： （2-6） 这里t是一个参数。也叫线性分式变换，因为当乘开之后得到这在每个变量上都是线性的，如果另一个固定的话；或者说在s和z上是双线性的。（1）在（2-11）中利用s=+j得到 （2-7）（2）整个左半平面映射到单位圆内，所以这是一个稳定的变换。（3）虚轴映射到单位圆是以一种一对一的方式进

6、行的，因此，在频域不存在混叠现象。在（2-9）中代入=0，得到（幅度为1）（2-8） 作为的函数解出得到 （2-9） 这表明是与非线性关联的，但是不存在混叠。因此，在（2-8）式中我们说是被预畴变到的。3 设计步骤3.1 设计流程图本次课程设计设计流程如下开始设计好的滤波器进行滤波处理验证并进行频谱分析语音信号的采集（wavread函数），画时域图快速傅里叶变换，并且画频谱图比较滤波前后语音信号的波形及频谱双线性变换法设计切比雪夫i型滤波器回放语音信号设定滤波器性能指标，通带截止频率fb=1700，阻带截止频率fc=1750，通带波纹ap=1，阻带波纹as=15结束图3-1 双线性变换法切比雪

7、夫i型滤波器对语音信号去噪流程图3.2 录制语音信号利用windows下的录音机，录制语音信号“今天，你吃了吗”，时间在1 s左右，要求为8000hz，8位单声道的音频格式。图3-2 window自带录音机图3-3 语音信号设置3.3 语音信号的频谱分析在matlab中编辑m函数，使用wavread函数读取采集的声音文件（.wav)它赋值给某一向量，再对其进行采样,然后使用plot语句画出相关的频谱图形（1）wavread函数调用格式：y,fs,nbits=wavread(file)（2）快速傅里叶变换算法fft计算dft的函数fft,其调用格式如下：x=fft(x,n)参数x为被变换的时域序

8、列向量，n是dft变换区间长度，当n大于x的长度时，fft函数自动在x后面补零。，当n小于xn的长度时，fft函数计算x的前n个元素，忽略其后面的元素。（3）声音采样文件读取的程序（文件名：liu.wav) ，如wavread(liu.wav)具体程序如下%语音信号的采集x,fs,bits=wavread(d:liu.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放 n=length(x); % 计算信号x的长度 fn=

9、2000; % 单频噪声频率，此参数可改 t=0:1/fs:(n-1)/fs; % 计算时间范围，样本数除以采样频率 x=x; y=x+0.05*sin(fn*2*pi*t); sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声 %语音信号的频谱分析 x=abs(fft(x); y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱 x=x(1:n/2); y=y(1:n/2); % 截取前半部分 deltaf=fs/n; % 计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围运行上述程序，得到图3.4 图3-

10、4 语音信号的时域图和fft频谱3.4 滤波器设计 设计指标：通带截止频率为1500hz，阻带截止频率为1550hz，通带波纹为1db，阻带波纹为15db，用脉冲（冲激）响应不变法设计的一个满足上述指标的切比雪夫i型滤波器 fb=1700;fc=1750; %定义通带和阻带截止频率 as=15;rp=1; %定义通带波纹和阻带衰减 wp=2*pi*fb/fs;ws=2*pi*fc/fs; %计算对应的数字频率 t=2; %定义采样间隔%双线性变换法设计% 第一步 预畸变 omegap=(2/t)*tan(wp /2);omegas=(2/t)*tan(ws /2); %截止频率预畸变% 第二步

11、 根据指标设计模拟滤波器 c,d=afd_chb1(omegap, omegas,rp,as);* 切比雪夫-1 滤波器阶次 = 13 b,a=bilinear(c,d,1/t); db,mag,pha,grd,w=freqz_m(b,a);运行上述程序，得到数字切比雪夫i型滤波器特性图如下图3-5 利用双线性变换设计的数字切比雪夫i型滤波器（w（单位：）3.4 信号滤波处理用双线性变换法设计好切比雪夫i型滤波器，那么我们就要对语音信号进行滤波处理。使用filter函数进行滤波，即y_fil = filter(b,a,y)，signal即是滤波后的语音信号，我们将滤波前后的时域波形进行对比，看

12、是否有处理好，并且我们对其进行快速傅里叶变换，即x=fft（y_fil），目的是对比滤波前后的频域频谱，进一步具体分析我们设计的滤波器是否达到我们的设计要求。 y_fil=filter(b,a,y); f_fil=abs(fft(y_fil); f_fil=f_fil(1:n/2);我们将滤波前后的语音信号进行时域对比，滤波前后频域频谱图如图3-6所示图3-6 原始语音信号及滤波前后的图形3.5 结果分析由图3-4可以看出，我们设计的阻带截止频率为0.525，此处的纵坐标差不多恰好在as处，所以设计的滤波器达到设计要求。由图3-5可以看出，滤波前后时域变化得不明显，这大概语音信号加入的噪声不大

13、，所以滤除噪声表现不明显；通过前后频域波形对比，可以很明显看到设计的滤波器对语音信号2000hz左右的噪声进行了滤除。另外通过sound函数听取原声和还原后的声音发现二者差别几乎难以识别。所以，双线性变换法设计的切比雪夫i型滤波器达到了设计要求。3.6 不同信源相同滤波方法的同学比较各种信源的特点图3-7为王晸同学所做的使用双线性变换法设计切比雪夫i型滤波器的音乐信号滤波去噪图3-7 原始音乐信号及滤波前后的图形图3-7与图3-5进行对比可以发现，幅度上难以看出差别，而图3-7中还原后音乐信号幅度谱与原始音乐信号幅度谱相差较大，因而相同滤波方法对音乐造成的失真较严重，分别运用sound函数可以

14、较明显听出二者差别，而语音信号则听起来差别不大。3.7 相同信源不同滤波方法的同学比较各种滤波方法性能借用李运达同学的双线性变换法设计的巴特沃斯滤波器模型得到图3-8 原始语音信号及滤波前后的图形（巴特沃斯滤波器）图3-8与图3-6进行对比，可见巴特沃斯滤波器在通带内幅度特性是单调下降的，如果阶次一定，在靠近截止频率处幅度下降很多，为了使通带内的衰减足够小，需要的阶次很高4 出现的问题及解决方法4.1 出现的问题（1）在作图时，matlab提示plot函数非法利用，向量要有相同的长度，因而无法作图。（2）滤波后得到的频率谱为空白，即得不到想要的滤除噪声的频率谱图。（3）得到的滤波后频率谱图中单

15、频噪声未曾滤去，声音信号有明显尖锐的单频啸叫声.（4）观察通带截止频率和阻带截止频率时，看得不够精准。4.2 解决方法（1）经过检查源程序发现，在程序中间非法插入fs=5000，与信号固有频率不同。（2）将通带截止频率该为1500hz，阻带截止频率该为1550hz，重新尝试，终于得出正确的图。（3）经过调试，更改合理的的参数，得到声音信号无明显噪声存在，用sound函数放出的还原声与原信号声相差不大。（4）将角频率的频谱图转化为频率的频谱图，这样在观察通带截止频率和阻带截止频率更加精确，误差更小。5 结束语本次课程设计总计两个星期，采用matlab进行数字信号处理和滤波器的设计。本次课程设计使

16、我加深了对课堂所掌握的概念理解，巩固并挖掘课堂上所学的各种知识，充分的意识到了我能做到些什么。采用matlab设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地。matlab 信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具, 在设计数字滤波器时, 善于应用matlab进行辅助设计, 能够大大提高设计效率。刚开始学习这门课程时,总感到数字信号处理的概念抽象,其中的分析方法与基本理论很难理解与掌握。但是此次课程设计,帮助我理解与掌握了课程中的基本概念、基本原理、基本分析方法以及综合应用所学知识解决实际问题的能力,深刻的证明了实践才能更快更

17、好的掌握一门繁琐的学科。在此我要特别感谢我们的指导老师，还有各位帮助我的同学，你们的指点教诲让我能在遇到困难及时的解决掉，无从下手时很快的找到突破口。在你们的帮助下，我将更加努力进步学习和完善自己。参考文献1 张圣勤. matlab7.0实用教程m. 北京：机械工业出版社，2006.32 维纳.k.英格尔，约翰.g.普罗克斯. 数字信号处理（matlab版）m. 西安：西安交通大学出版社，2008.13 陈金鹰.dsp技术及应用m. 北京：机械工业出版社.2004.64 孙宗瀛. dsp原理设计与应用m. 北京：清华大学出版社，20025 薛年喜主编.matlab 在数字信号处理中的应用m北京

18、：清华大学出版社 附录1：源程序清单% 程序功能：采集一段语音信号，绘制波形并观察其频谱，设计一个切比雪夫i型滤波器，并对语音信号进行滤波去噪% 程序作者：刘新武% 最后修改日期：2010-03-04%语音信号的采集x,fs,bits=wavread(d:liu.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放 n=length(x); % 计算信号x的长度 fn=2000; % 单频噪声频率，此参数可改 t=0:1/f

19、s:(n-1)/fs; % 计算时间范围，样本数除以采样频率 x=x; y=x+0.05*sin(fn*2*pi*t); sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声 %语音信号的频谱分析 x=abs(fft(x); y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱 x=x(1:n/2); y=y(1:n/2); % 截取前半部分 deltaf=fs/n; % 计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围 subplot(2,2,1);plot(t,x); title(原始语音信号);xlab

20、el(时间 );ylabel(音量 );%做原始语音信号幅度 subplot(2,2,2);plot(x);title(语音信号幅度谱图); subplot(2,2,3);plot(t,y); title(加入单频干扰后语音信号);xlabel(时间 );ylabel(音量 );%做加入单频干扰后语音信号 subplot(2,2,4);plot(y);title(加入单频干扰后语音信号幅度谱图); fb=1700;fc=1750; %定义通带和阻带截止频率 as=15;rp=1; %定义通带波纹和阻带衰减 wp=2*pi*fb/fs;ws=2*pi*fc/fs; %计算对应的数字频率 t=2;

21、 %定义采样间隔%双线性变换法设计% 第一步 预畸变 omegap=(2/t)*tan(wp /2);omegas=(2/t)*tan(ws /2); %截止频率预畸变% 第二步 根据指标设计模拟滤波器 c,d=afd_chb1(omegap, omegas,rp,as);* 切比雪夫-1 滤波器阶次 = 21 b,a=bilinear(c,d,1/t); db,mag,pha,grd,w=freqz_m(b,a); subplot(2,2,1);plot(w,db) ;title(幅度响应);xlabel(w/pi);ylabel(|h|); subplot(2,2,2);plot(w,ma

22、g) ;title(幅度（db）);xlabel(w/pi);ylabel(db); subplot(2,2,3);plot(w,grd) ;title(相位响应);xlabel(rad);ylabel(h（n）); subplot(2,2,4);plot(t,pha) ;title(脉冲响应);xlabel(rad);ylabel(h（t）);y_fil=filter(b,a,y); %进行滤波处理 f_fil=abs(fft(y_fil); f_fil=f_fil(1:n/2); subplot(3,2,1); plot(t,x);title(原始语音信号);xlabel(时间);ylab

23、el(音量); % 画出原始语音信号时域图 subplot(3,2,2); plot(f,x); title(原始语音信号幅度谱);xlabel(频率f);ylabel(幅度);% 画出原始语音信号幅度谱 subplot(3,2,3); plot(t,y);title(加入噪声后语音信号); xlabel(时间);ylabel(音量);% 画出加入噪声后语音信号时域图 subplot(3,2,4); plot(f,y);title(加入噪声后语音信号幅度谱);xlabel(频率f);ylabel(幅度);% 画出加入噪声后语音信号幅度谱 subplot(3,2,5); plot(t,y_fil

24、);title(滤去噪声后语音信号);xlabel(时间);ylabel(音量);% 画出滤去噪声后语音信号时域图 subplot(3,2,6); plot(f,f_fil); title(滤去噪声后语音信号幅度谱) ;xlabel(频率f);ylabel(幅度);% 画出滤去噪声后语音信号幅度谱附录2：自编程函数清单function b,a = afd_chb1(wp,ws,rp,as);% 切比雪夫-1型模拟低通滤波器设计% -% b,a = afd_chb1(wp,ws,rp,as);% b = ha(s) 分子的系数% a = ha(s) 分母的系数% wp = 以弧度/秒为单位的通带边缘频率; wp 0% ws = 以弧度/秒为单位的阻带边缘频率; ws wp 0% rp = 通带中的振幅波动的+db数; (rp 0)% as = 阻带衰减的+db数; (as 0)%if wp = 0 error(通带边缘必须大于 0)endif ws = wp error(阻带边缘必须大于通带