语音信号滤波去噪——使用双线性变换法设....doc

《语音信号滤波去噪——使用双线性变换法设....doc》由会员分享，可在线阅读，更多相关《语音信号滤波去噪——使用双线性变换法设....doc（19页珍藏版）》请在一课资料网上搜索。

1、 语音信号滤波去噪使用双线性变换法设计的切比雪夫I型滤波器 第 页 共19页语音信号滤波去噪使用双线性变换法设计的并联型切比雪夫I型滤波器学生姓名：欧晓燕 指导老师：吴志敏摘 要 本课程设计是采用双线性变换法设计的切比雪夫I型滤波器对双音频信号滤波去噪。在网上下载一段双音频信号，在MATLAB集成环境下，首先用wavread函数求出双音频信号的相关参数，对双音频信号进行读取和加噪；然后再给定相应技术指标，设计一个满足指标的切比雪夫I型滤波器，对该双音频信号进行滤波去噪处理，并绘制对比图，比较滤波前后的波形和频谱并进行分析；最后通过回放双音频信号，对比滤波前后的信号变换。本课程设计成功的对双音频

2、信号进行滤波去噪，初步完成了设计指标。关键词 双音频信号； 滤波设计； MATLAB；切比雪夫I型滤波器1 引 言用麦克风采集一段8000Hz，8k的双音频信号，绘制波形并观察其频谱，给定通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指标的切比雪夫I型IIR滤波器，对该双音频信号进行滤波去噪处理。1.1 课程设计目的数字信号处理课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实

3、验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2 课程设计的要求（1）学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法；（2）滤波器指标必须符合工程实际，根据模拟滤波器的性能指标，确定数字滤波器指标；（3）采用双线性变换法，设计满足上述性能指标要求的ChebyshevI型数字低通滤波器；（4）设计完后应检查其频率响应曲线是否满足指标；（5）处理结果和分析结论应该一致，而且应符合理论；（6）独立完成课程设计并按要求编写课程设计报告书；

4、1.3 设计平台本次课程设计是在MATLAB软件平台上进行的。MATLAB是矩阵实验室（MATRIX LABORATORY）的简称，是美国MATHWORKS公司推出的具有强大数值分析、矩阵运算、图形绘制和数据处理等功能的软件，现已广泛应用到教学、科研、工程设计等领域2。随着MATLAB软件信号处理工具箱的推出，MATLAB已成为信息处理，特别是数字信号处理DSP应用中分析和设计的主要工具。就MATLAB信号处理中的滤波器设计而言，在很大程度上能快速有效地实现滤波器的分析、设计及仿真，大大节约了设计时间，相对传统设计而言，简化了滤波器设计的难度。2 设计原理用麦克风采集一段双音频信号，绘制波形并

5、观察其频谱，给定相应技术指标，用双线性变换法设计的切比雪夫I型IIR滤波器，对该双音频信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 IIR滤波器从离散时间来看，若系统的单位抽样(冲激)响应延伸到无穷长，称之为“无限长单位冲激响应系统”，简称为IIR系统。无限长单位冲激响应（IIR）滤波器有以下几个特点:(1) 系统的单位冲激响应h(n)是无限长；(2) 系统函数H(z)在有限z平面（0）；(3) 结构上存在着输出到输入的反馈，也就是结构上是递归型的。IIR滤波器采用递归型结构，即结构上带有反馈环路。同一种系统函数H（z）可以有多种不同的结构，基本网络结构有直接型、直接型、级

6、联型、并联型四种，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些典型的滤波器各有特点。有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。2.2 切比雪夫I型滤器切比雪夫滤波器（又译车比雪夫滤波器）是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”

7、。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。这种滤波器来自切比雪夫多项式，因此得名，用以纪念俄罗斯数学家巴夫尼提列波维其切比雪夫 巴特沃兹滤波器在通带内幅度特性是单调下降的，如果阶次一定，则在靠近截止 处，幅度下降很多，或者说，为了使通带内的衰减足够小，需要的阶次N很高，为了克服这一缺点，采用切比雪夫多项式来逼近所希望的 。切比雪夫滤波器的 在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。切比雪夫滤波器的振幅平方函数为 （2-1） 式中

8、 c有效通带截止频率 与通带波纹有关的参量， 大,波纹大 0 1时, |x|, VN(x) 切比雪夫滤波器的振幅平方特性如图所示,通带内的变化范围为 1(max) (min)时，|x|1,随 ， 0 (迅速趋于零)当 =0时， （2-3） N为偶数，cos2( )=1，得到min， （2-4）N为奇数，cos2( ，得到max， （2-5）切比雪夫滤波器的振幅平方特性如图2.1所示。图2.1 切比雪夫滤波器的振幅平方特性2.3 双线性变换法双线性变换法是使数字信号滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了客服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面

9、的一条横带里（宽度为,即从到),其次再通过上面讨论过的标准变换关系将此横带变换到这个z平面上去，这样就使s平面与z平面式一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。将s平面整个轴压缩变换到s1平面轴上的到一段，可以采用以下变换关系： （2-6） 这样，变换到，变到可将（6）式写成 （2-7）解析延拓到整个s平面和s1平面，令，则得 （2-8）再将s1平面通过以下标准变换关系映射到z平面： （2-9）从而得到s平面和z平面的单值映射的关系为 （2-10） （2-11）一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应的关系，可以引入待定常数c，使（6）式和（7）式变换

10、成 （2-12） （2-13）仍将代入（13）式，可得 （2-14） （2-15） （14）式和（15）式是s平面与z平面之间的单值映射关系，这种变换就称为双线性变+换。3 设计步骤3.1设计流程图语音信号滤波去噪使用双线性变换法设计的切比雪夫I型滤波器的设计流程如图3.1所示：开始结束回放音频信号比较滤波前后双音频信号的波形及频谱设计好的滤波器进行滤波处理验证并进行频谱分析语音信号的采集（wavread函数），画时域图快速傅里叶变换，并且画频谱图设定滤波器性能指标，通带截止频率fb=2000，阻带截止频率fc=2150，通带波纹Rp=1，阻带波纹As=35双线性变换法设计切比雪夫I型滤波器图

11、3.1 双线性变换法切比雪夫I型滤波器对双音频信号去噪流程图3.2语音信号的采集点击windows系统桌面的“开始”按钮，点击开始菜单栏里的“附件”，选择“录音机”选项，点击录音机“文件”选项，进入“声音选定”设置，把属性一栏设置成“8000Hz，8位，单声道，7KB/秒”(见图3.2)。点击确定，然后开始双音频信号的采集，采集时间为23秒左右为最佳。采集的声音文件以“.wav”格式存储（见图3.3）。图3.2 采集声音的参数设置图3.3 采集声音3.3语音信号的频谱分析在MATLAB中编辑m函数，使用wavread函数读取采集的声音文件（.wav)将它赋值给某一向量，再对其进行采样,然后使用

12、plot语句画出相关的频谱图形。（1）Wavread函数调用格式：y,fs,nbits=wavread(file)功能说明：采样值放在向量x中，fs表示采样频率（Hz），bits表示采样位数。（2）快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：Xk=fft(x,n)参数x为被变换的时域序列向量，N是DFT变换区间长度，当n大于x的长度时，fft函数自动在x后面补零。，当n小于xn的长度时，fft函数计算x的前n个元素，忽略其后面的元素。在本次课程设计中，我们利用fft函数对双音频信号进行快速傅里叶变换，就可以得到信号的频谱特性。（3）声音采样文件读取的程序（文件名：xiaom

13、iao.wav)；双音频信号的提取：x,fs,bits=wavread(D:MATLAB7workxiaomiao.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2200; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); y=x+sin(fn*2*pi*t);sound(y,fs,bits); %

14、 应该可以明显听出有尖锐的单频啸叫声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.4 原始双音频信号和加噪后信号的时域和频谱图3.4滤波器设计 设计指标：通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指

15、标的切比雪夫I型滤波器双线性变换法设计切比雪夫I型滤波器fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=35; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi; ws=fc/fs*2*pi; %计算对应的数字频率T=1;Fs = 1/T; %定义采样间隔Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2); %截止频率预畸变cs,ds = afd_chb1(Omegap,Omegas,Rp,As); %选择滤波器最小阶数b,a = bilinear(cs,ds,Fs);C,B,A = dir2cas(b,a) %双线性

16、变换法实现模拟滤波器到数字滤波器的转换db,mag,pha,grd,w=freqz_m(b,a); %绘制数字滤波器频率响应幅度图delta=1,zeros(1,99);ha=filter(b,a,delta);得到切比雪夫滤波器的幅度和相位谱如图3.5所示。图3.5 切比雪夫滤波器的幅度和相位谱3.5双音频的滤波程序及滤波效果图滤波程序：y_fil=filter(b,a,y); % 用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半画出滤波前后双音频信号的时域和频谱图如图3.6所示。图3.6 滤波前后双音频信号的

17、时域和频谱图3.6结果分析在MATLAB中，经sound函数，对经过切比雪夫I型滤波器之后的信号进行回放，可以听出滤波之后的信号比原始信号更清晰一些，清除了环境噪音。通过以下语句来进行语音信号回放比较：sound（x，fs，bits） %播放原始双音频信号sound（y，fs，bits）%播放加噪后的双音频信号sound（y_fil，fs，bits） %播放经过滤波处理后的双音频信号所得结果证明了切比雪夫I型滤波器去噪设计成功。4 出现的问题及解决方法 在这次的课程设计中，由于理论知识的不踏实以及其他各种原因，我们遇到了不少问题。（1）在进行双音频信号提取时，进过多次录取才得到理想的双音频信号

18、，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的双音频信号的时域图和频谱图（2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的双音频信号wav文件放到Matlab的work文件夹中。（3）还要在滤波器性能曲线的As处画一根竖线，这样更方便看出结果。（4）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。5 课程设计心得体会采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是MATLAB

19、工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。MATLAB 信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具, 在设计数字滤波器时, 善于应用MATLAB进行辅助设计, 能够大大提高设计效率。两周的课程设计让我对MATLAB软件的使用更加的熟练，对切比雪夫滤波器的滤波原理有了更深刻的认识。6 参考文献1 程佩青数字信号处理教程M北京：清华大学出版社，20022 薛年喜主编MATLAB 在数字信号处理中的应用M北京：清华大学出版社，20023 维纳K恩格尔，约翰G普罗克斯.数字

20、信号处理M.西安交通大学出版社，20024 董长虹等. MATLAB信号处理与应用M.北京：国防工业出版社，20055 美 M.H.海因斯 著，张建华等 译.数字信号处理M.北京：科学出版社，20026 张葛祥，李 娜. MATLAB仿真技术与应用M.北京：清华大学出版社，2007附录一: 源程序设计清单原始双音频信号的采集及加噪：x,fs,bits=wavread(D:MATLAB7workxiaomiao.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); %

21、 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2200; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); y=x+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); % 截取前半部分Warning: Integer operands are required for colon

22、 operator when used as index.Warning: Integer operands are required for colon operator when used as index.deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围subplot(2,2,1);plot(t,x);xlabel(时间(单位:s);ylabel(幅度);title(原始双音频信号);gridsubplot(2,2,2);plot(f,X);xlabel(频率(单位:Hz);ylabel(幅度谱);title(原始双音频

23、信号幅度谱图); axis(0,6*103,0,6000);gridsubplot(2,2,3);plot(t,y);xlabel(时间(单位:s);ylabel(幅度);title(加噪后的双音频信号);gridsubplot(2,2,4);plot(f,Y);xlabel(频率(单位:Hz);ylabel(幅度谱);title(加噪后的双音频信号幅度谱图); axis(0,6*103,0,6000);grid使用双线性变换法设计切比雪夫I型滤波器：fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=35; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;

24、ws=fc/fs*2*pi; %计算对应的数字频率T=1;Fs = 1/T; %定义采样间隔Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2); %截止频率预畸变cs,ds = afd_chb1(Omegap,Omegas,Rp,As); %选择滤波器最小阶数Warning: Function call afd_chb1 invokes inexact match d:MATLAB7workAFD_CHB1.M.* 切比雪夫-1 滤波器阶次 = 14 Warning: Function call u_chb1ap invokes inexact matc

25、h d:MATLAB7workU_CHB1AP.M. In AFD_CHB1 at 29b,a = bilinear(cs,ds,Fs);C,B,A = dir2cas(b,a) %双线性变换法实现模拟滤波器到数字滤波器的转换Warning: Function call dir2cas invokes inexact match d:MATLAB7workDIR2CAS.M.C =3.3307e-016B = 1.0000 -0.0321 0.4021 1.0000 -0.2091 0.0238 1.0000 -0.2206 5.8342 1.0000 -0.3400 0.1487 1.000

26、0 -0.5438 0.7109 1.0000 -2.1153 8.5100 1.0000 -12.5391 364.1406A = 1.0000 -1.9050 0.9776 1.0000 -1.9116 0.9540 1.0000 -1.9147 0.9950 1.0000 -1.9185 0.9754 1.0000 -1.9315 0.9425 1.0000 -1.9354 0.9570 1.0000 -1.9430 0.9444db,mag,pha,grd,w=freqz_m(b,a); %绘制数字滤波器频率响应幅度图Warning: Function call freqz_m inv

27、okes inexact match d:MATLAB7workFREQZ_M.M. delta=1,zeros(1,99);ha=filter(b,a,delta);figureSubplot(221);plot(w/pi,db); title(切比雪夫I 滤波器幅度(dB);xlabel(w);ylabel(dB);axis(0,1,-150,10);gridSubplot(222);plot(w/pi,mag);title(切比雪夫I滤波器幅度响应);xlabel(w);ylabel(幅值 |H|); axis(0,1,0,1);gridSubplot(223);plot(w/pi,ph

28、a);title(切比雪夫I 滤波器相位响应); xlabel(w);ylabel(相位（单位： ）);axis(0,1,-4,4);grid对加噪信号进行滤波处理并画出加噪前后的信号时域和频谱：y_fil=filter(b,a,y); % 用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半Warning: Integer operands are required for colon operator when used as index.figuresubplot(3,2,1);plot(t,x);xlabel

29、(时间(单位:s);ylabel(幅度);title(原始双音频信号);gridsubplot(3,2,2);plot(f,X);xlabel(频率(单位:Hz);ylabel(幅度谱);title(原始双音频信号幅度谱图);gridaxis(0,6*103,0,6000)subplot(3,2,3);plot(t,y);xlabel(时间(单位:s);ylabel(幅度);title(加噪后的双音频信号);gridsubplot(3,2,4);plot(f,Y);xlabel(频率(单位:Hz);ylabel(幅度谱);title(加噪后的双音频信号幅度谱图);gridaxis(0,6*10

30、3,0,6000)subplot(3,2,5);plot(t,y_fil);xlabel(时间(单位:s);ylabel(幅度);title(滤波后的双音频信号);gridaxis(0,6,-1,1)subplot(3,2,6);plot(f,Y_fil);xlabel(频率(单位:Hz);ylabel(幅度谱);title(滤波后的双音频信号幅度谱图);gridaxis(0,6*103,0,6000)附录二：afd_chb1函数function b,a = afd_chb1(Wp,Ws,Rp,As);% 切比雪夫-1型模拟低通滤波器设计% -% b,a = afd_chb1(Wp,Ws,Rp

31、,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(阻带边缘必须大于通带边缘)endif (Rp = 0) | (As 0) error(通带波动或阻带衰减必须大于0)endep = sqrt(10(Rp/10)-1);A = 10(As/20);Omeg

32、aC = Wp;OmegaR = Ws/Wp;g = sqrt(A*A-1)/ep;N = ceil(log10(g+sqrt(g*g-1)/log10(OmegaR+sqrt(OmegaR*OmegaR-1);fprintf(n* 切比雪夫-1 滤波器阶次 = %2.0f n,N)b,a=u_chb1ap(N,Rp,OmegaC);附录三：freqz.m.m函数function db,mag,pha,grd,w = freqz_m(b,a);% freqz 子程序的改进版本% -% db,mag,pha,grd,w = freqz_m(b,a);% db = 0 到pi弧度区间内的相对振幅(

33、db)% mag = 0 到pi弧度区间内的绝对振幅% pha = 0 到pi弧度区间内的相位响应% grd = 0 到pi弧度区间内的群迟延% w = 0 到pi弧度区间内的501个频率样本向量% b = Ha(z)的分子多项式系数(对FIR b=h)% a = Ha(z)的分母多项式系数(对 FIR: a=1)%H,w = freqz(b,a,1000,whole); H = (H(1:1:501); w = (w(1:1:501); mag = abs(H); db = 20*log10(mag+eps)/max(mag); pha = angle(H);% pha = unwrap(angle(H); grd = grpdelay(b,a,w);% grd = diff(pha);% grd = grd(1) grd;% grd = 0 grd(1:1:500); grd; grd(2:1:501) 0;% grd = median(grd)*500/pi;