数字信号课程设计：语音信号的采集、分析与处理

1、长沙学院长 沙 学 院数字信号处理课程设计说明书题 目数字信号处理课程设计系（部）电子与通信工程系专业（班级）10级通信2班姓 名周斌学 号2010043205指导老师李广柱、刘光灿、陈威兵黄飞江、张刚林、冯 璐起止日期2013.4.15-2013.4.26- 27 -电子与通信工程系目 录1.绪论- 7 -2.设计作品名称- 7 -3.课程设计内容- 7 -3.1 设计思想- 7 -3.2设计原理- 8 -3.2.1 无限脉冲响应数字滤波器设计- 8 -3.2.1.1 巴特沃斯滤波器设计：- 8 -3.1.1.2 巴特沃斯滤波器的特性：- 8 -3.1.1.3 双线性变换法：- 10 -3.

2、2.2 有限脉冲响应数字滤波器设计- 12 -4. 课程设计设计步骤及结果分析- 15 -4.1 语音信号的采集及其时域、频域分析- 15 -4.1.1程序- 15 -4.1.2 图像- 15 -4.2 设计数字滤波器和画出其频率相映- 16 -4.2.1 无限低通数字滤波器- 16 -4.2.1.1 程序- 16 -4.2.1.2 图像- 17 -4.2.2 无限高通数字滤波器- 17 -4.2.2.1 程序- 17 -4.2.2.2 图像- 18 -4.2.3 有限低通数字滤波器- 19 -4.2.3.1程序- 19 -4.2.3.2 图像- 19 -4.2.4 有限高通数字滤波器- 20

3、 -4.2.4.1 程序- 20 -4.2.4.2 图像- 20 -4.3 用滤波器对信号进行滤波并比较前后波形- 21 -4.3.1 经过低通iir- 21 -4.3.2 经过高通iir- 22 -4.3.3 经过低通fir- 23 -4.3.4 经过高通fir- 24 -4.4 编制实现上述任务的相应的总程序- 24 -4.4.1 信号经过低通iir- 24 -4.4.2 信号经过高通iir- 25 -4.4.3 信号经过低通fir- 27 -4.4.4 信号经过高通 fir- 28 -4.5 撰写5000-8000字课程设计报告- 29 -5. 总结- 29 -6. 存在的建议及不足-

4、30 -7. 参考文献- 30 -摘 要 对一段语音信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。 通过数字信号处理的课程设计，使学生对信号的采集，处理，传输，显示，存储和分析等有一个系统的掌握和理解。巩固和运用数字信号处理课程中的理论知识和实验技能，掌握最基本的数字信号处理的理论和方法，培养学生发现问题，分析问题和解决问题的能力。为以后的工作

5、和学习打下基础。本课题主要基于matlab平台，实现信号分析与滤波器的设计。并改变滤波器参数或特性(低通、高通)，实现不同的滤波要求。自己录播一段语音信号，对其进行采样，进行频谱分析，设计有限和无限两种不同的数字滤波器对信号进行滤波处理，分析所设计滤波器(画出了频率特性曲线)，并对信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化，分析的结果表明初步实现了设计目标。关键词：无限脉冲响应数字滤波器、有限脉冲响应数字滤波器、巴特沃斯滤波器、窗函数法、双线性法 1.绪论 对一段语音信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函

6、数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。 通过数字信号处理的课程设计，使学生对信号的采集，处理，传输，显示，存储和分析等有一个系统的掌握和理解。巩固和运用数字信号处理课程中的理论知识和实验技能，掌握最基本的数字信号处理的理论和方法，培养学生发现问题，分析问题和解决问题的能力。为以后的工作和学习打下基础。本课题主要基于matlab平台，实现信号分析与滤波器的设计。并改变滤波器参数或特性(低通、高通)，实现不同的滤波要求。

7、自己录播一段语音信号，对其进行采样，进行频谱分析，设计有限和无限两种不同的数字滤波器对信号进行滤波处理，分析所设计滤波器(画出了频率特性曲线)，并对信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化，分析的结果表明初步实现了设计目标。2.设计作品名称语音信号的采集、分析与处理3.课程设计内容3.1 设计思想采样一段语音信号；画出语音信号的时域波形和频谱图；给定滤波器的性能指标，设计数字滤波器，并画出滤波器的频率响应；然后用设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱。关键技术：频谱图的理解；设计数字滤波器；数字滤波的方法；解决思路：对语

8、音号进行快速傅里叶变换，得到信号的频谱特性；在matlab环境中可以利用函数fir1设计fir滤波器，可以利用函数butter，cheby1和ellip设计iir滤波器；利用matlab中的函数freqz画出各滤波器的频率响应。3.2设计原理3.2.1 无限脉冲响应数字滤波器设计3.2.1.1 巴特沃斯滤波器设计：巴特沃斯滤波器的特点是同频带的频率响应去想最平滑，如下图所示图1 一级巴特沃斯低通滤波器的波得图图2 一级至五级巴特沃斯低通滤波器3.1.1.2 巴特沃斯滤波器的特性：巴特沃斯滤波器的特点是同频带内的频率响应曲线最为平坦，没有起伏，而在组频带则逐渐下降为零。在振幅的对数对角频率的波特

9、图上，从某一边界见频率开始，振幅随着角频率的增加而逐渐减少，趋向于负无穷大。一阶巴特沃斯滤波器的衰减率为每倍频20分贝，二阶巴特沃斯滤波器的衰减率为每倍频12分贝，三阶的衰减率为每分贝18分贝，如此类推，巴特沃斯滤波器的振幅对角频率单调下降，并且滤波器的结束越高，在组频带振幅衰减速度越快，其他滤波器高阶的振幅对角频率图和低阶数的振幅对角频率有不同的形状。 上述函数的特点是等距离分布在半径为的圆上。因此，极点用下式表示为的表示式：为了使设计公式和图表统一，将频率归一化。巴特沃斯滤波器采用3db截止频率归一化，归一化后的系统函数为令,称为归一化频率，称为归一化复变量，这样巴特沃斯滤波器的归一化低通

10、原型系统函数为式中，为归一化极点，用下式表示: 然后将巴特沃斯滤波器转换成高通滤波器，转换公式如下：在平面与平面虚轴上的频率关系为：其中，为希望设计的高通滤波器的通带边界频率。3.1.1.3 双线性变换法：双线性变换法是将平面压缩变换到某一中介平面的一条横带里，再通过标准变换关系将此带变换到整个z平面上去，这样就使平面与平面之间建立一一对应的单值关系，消除了多值变换性。为了将平面的轴压缩到平面的轴上的到一段上，可以通过以下的正切变换来实现：这样当由经变化到时，由经过0变化到，也映射到了整个 轴。将这个关系延拓到整个平面和平面，则可以得到再将平面通过标准变换关系映射到平面，即令得到同样对z求解，

11、得到这样的变换叫做双线性变换。为了验证这种映射具有s平面的虚轴映射到z平面单位圆上的特性，考虑 ,，得除了使s平面的虚轴映射到单位圆上之外，s平面的左半部分映射到单位圆的内部，s平面的右半部分映射到单位圆的外部。如图所示图3 双线性变化映射关系示意图观察式子，发现的实部为负时，因子的幅度小于1，相当于单位圆的内部。反之，当的实部为负时，该比值的幅度大于1，相当于单位圆的外部。这样就可以看出使用双线性变换可从稳定的模拟滤波器得到稳定的数字滤波器。双线性变换法还避免了使用脉冲响应不变法所遇到的混叠问题，因为它把平面的这个虚轴映射到平面的单位圆上。然而，付出的代价是在频率轴上引入了失真。因此，只有当

12、能容忍或补偿这种失真时，使用双线性变换法设计数字滤波器的方法才是实用的。仅在零频率附近时与之间的频率变换关系接近于线性关系，所产生的数字滤波器的幅频响应相对于原模拟滤波器的幅频响应有畸变。对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸变来加以校正，也就是将临界频率事先加以畸变，然后经变换后正好映射到所需要的频率上。通过的关系变换成一组模拟频率。图4 双线性变化法的频率关系为了克服冲击响应不变法产生的频率混叠现象，我们需要使平面与平面建立一一对应的单值关系，即求出，然后将其代入就可以求得，即 3.2.

13、2 有限脉冲响应数字滤波器设计如果所希望的滤波器的理想的频率响应函数为，则其对应的单位脉冲响应为 （4.1）窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，而且是非因果的，所以用窗函数将截断，并进行加权处理，得到：（4.2）就作为实际设计的fir数字滤波器的单位脉冲响应序列，其频率响应函数为（4.3）式中，n为所选窗函数的长度。我们知道，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度n的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度n 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。表1 各种窗函数的基本

14、参数窗函数旁瓣峰值幅度/db过渡带宽阻带最小衰减/db矩形窗-134/n-12三角形窗-258/n-25汉宁窗-318/n-44哈明窗-418/n-53不莱克曼窗-5712/n-74凯塞窗(=7.865)-5710/n-80这样选定窗函数类型和长度n之后，求出单位脉冲响应，并按照式（4.3）求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用fft计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度n ，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足根据上式中的正、负号和长度n的奇偶性又将线性相位fir滤波器分成四类。要根据所

15、设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。主程序框图如图4.1所示。其中幅度特性要求用db表示。开始读入窗口长度n计算hd(n)调用窗函数子程序求w(n)调用子程序（函数）计算h(k)=dfth(n)调用绘图子程序（函数）绘制h(k)幅度相位曲线结束图5 主程序框图计算h(n)= hd(n) w(n)设画图时，用打印幅度特性。第k点对应的频率。为使曲线包络更接近的幅度特性曲线，dft变换区间要选大些。例如窗口长度n=33时，可通过在末尾补零的方法，使长度变为64，再进行64点dft，则可以得到更精确的幅度衰减特性曲线。 4. 课程设计设计步

16、骤及结果分析4.1 语音信号的采集及其时域、频域分析对一段语音信号，时间在2s内。然后在matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解掌握采样频率、采样位数等概念。要求首先画出采样后语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。4.1.1程序x1,fs,bits=wavread(c:matlab7workcricket.wav);%原信号sound(x1,fs,bits); %调用sound函数播放原始语音figure(1); plot(x1);%做原始语音信号

17、的时域图形 title(原始语音信号时域图形);xlabel(时间 t); ylabel(音量 n);figure(2);y1=fft(x1);%做length(x1)点的fft y1=fftshift(y1); %平移是频率中心为0 derta_fs = fs/length(x1); %设置频谱的间隔分辨率 plot(-fs/2:derta_fs: fs/2-derta_fs,abs(y1);%画出原始语音信号的频谱图title(原始语音信号的频谱图); grid on;4.1.2 图像4.2 设计数字滤波器和画出其频率相映给出各滤波器的性能指标：(1) 低通滤波器性能指标 fb1000hz

18、，fc1200hz，as100db，ap1db。(2) 高通滤波器性能指标 fc4800hz，fb5000hz，as100db，ap1db。要求用窗函数法和双线性变换法设计上面要求的2种滤波器。在matlab中，可以利用函数fir1设计fir滤波器，可以利用函数butter设计iir滤波器；利用matlab中的函数freqz画出各滤波器的频率响应。4.2.1 无限低通数字滤波器4.2.1.1 程序wp=1000*2*pi;ws=2400*2*pi;rp=1;as=100;n,wc=buttord(wp,ws,rp,as,s);z,p,k=buttap(n);bap,aap=zp2tf(z,p,

19、k);bl,al=lp2lp(bap,aap,wc);b,a=bilinear(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(iir低通滤波器频率响应);4.2.1.2 图像4.2.2 无限高通数字滤波器4.2.2.1 程序fs=22050; %巴特沃斯高通滤波器wp=4800*2*pi;ws=10000*2*pi;rp=1;as=100;n,wc=buttord(wp,ws,rp,as,s);z,p,k=buttap(n);bap,aap=zp2tf(z,p,k);bl,al=lp2hp(bap,aap,wc);

20、b,a=bilinear(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(iir高通滤波器频率响应);4.2.2.2 图像4.2.3 有限低通数字滤波器4.2.3.1程序n=length(x1)-1;t=0:1/fs:n/fs;fp=1000;fc=1500;as=100;ap=1;fs=22050;%滤波器指标wc=2*pi*fc/fs; wp=2*pi*fp/fs;bt=wc-wp; %计算过渡带宽alph=0.112*(as-8.7); %alph用来控制窗的形状n=ceil(as-8)/2.285/bt);

21、%凯塞窗滤波器阶数n,(h(n)的长度为n+1)wn= kaiser(n+1,alph); ws=(wp+wc)/2/pi; %计算通带截止频率（关于pi归一化）hn=fir1(n,ws,wn);figure(2) %画出滤波器频谱图freqz(hn,1)4.2.3.2 图像4.2.4 有限高通数字滤波器4.2.4.1 程序n=length(x1)-1;t=0:1/fs:n/fs;fc=4800;fp=5000;as=100;ap=1;fs=22050; %采样频率必须变大（n-1）ws=2*pi*fc/fs; wp=2*pi*fp/fs;bt=wp-ws;n=ceil(8*pi/bt);wc

22、=(wp+ws)/2/pi;b,a=fir1(n,wc,high);figure(2);freqz(b,a);4.2.4.2 图像4.3 用滤波器对信号进行滤波并比较前后波形要求用自己设计的各种滤波器分别对采集的信号进行滤波，在matlab中，fir滤波器利用函数fftfilt对信号进行滤波，iir滤波器利用函数filter对信号进行滤波。4.3.1 经过低通iir4.3.2 经过高通iir4.3.3 经过低通fir4.3.4 经过高通fir4.4 编制实现上述任务的相应的总程序4.4.1 信号经过低通iir%无限脉冲响应数字低通滤波器（低通iir）clc;clear all;x1,fs,bi

23、ts=wavread(c:matlab7workcricket.wav);datlen=length(x1);sound(x1,fs,bits);y1=fftshift(fft(x1);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y1);title(原始语音信号频谱)wp=1000*2*pi;ws=2400*2*pi;rp=1;as=100;n,wc=buttord(wp,ws,rp,as,s);z,p,k=buttap(n);b

24、ap,aap=zp2tf(z,p,k);bl,al=lp2lp(bap,aap,wc);b,a=bilinear(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(iir低通滤波器频率响应);x=filter(b,a,x1); %将信号经过滤波器sound(x,fs,bits);figure(3) %画出滤波后的信号subplot(2,1,1)plot(x);title(过滤后语音信号时域波形);subplot(2,1,2)y3=fftshift(fft(x)plot(-datlen/2:datlen/2-1)/dat

25、len*fs,abs(y3);title(过滤后语音信号频域波形);4.4.2 信号经过高通iir%无限脉冲响应数字高通滤波器（高通iir）clear all;x1,fs,bits=wavread(c:matlab7workcricket.wav);datlen=length(x1);sound(x1,fs,bits);y1=fftshift(fft(x1);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y1);title(原始语音

26、信号频谱)fs=22050; %巴特沃斯高通滤波器wp=4800*2*pi;ws=10000*2*pi;rp=1;as=100;n,wc=buttord(wp,ws,rp,as,s);z,p,k=buttap(n);bap,aap=zp2tf(z,p,k);bl,al=lp2hp(bap,aap,wc);b,a=bilinear(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(iir高通滤波器频率响应);x=filter(b,a,x1); %将信号经过滤波器sound(x,fs,bits);figure(3) %画出

27、滤波后的信号subplot(2,1,1)plot(x);title(过滤后语音信号时域波形);subplot(2,1,2)y3=fftshift(fft(x)plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y3);title(过滤后语音信号频域波形);4.4.3 信号经过低通fir%信号经过用窗函数法设计低通firclc;clear all;x1,fs,bits=wavread(c:matlab7workcricket.wav);datlen=length(x1);sound(x1,fs,bits);y1=fftshift(fft(x1);figure(1)s

28、ubplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y1);title(原始语音信号频谱)n=length(x1)-1;t=0:1/fs:n/fs;fp=1000;fc=1500;as=100;ap=1;fs=22050;%滤波器指标wc=2*pi*fc/fs; wp=2*pi*fp/fs;bt=wc-wp; %计算过渡带宽alph=0.112*(as-8.7); %alph用来控制窗的形状n=ceil(as-8)/2.285/bt); %凯塞窗滤波器阶数n

29、,(h(n)的长度为n+1)wn= kaiser(n+1,alph); ws=(wp+wc)/2/pi; %计算通带截止频率（关于pi归一化）hn=fir1(n,ws,wn);figure(2) %画出滤波器频谱图freqz(hn,1)x=conv(hn,x1); %将信号经过滤波器sound(x,fs,bits);figure(3) %画出滤波后的信号subplot(2,1,1)plot(x);title(过滤后语音信号时域波形);subplot(2,1,2)y=fftshift(fft(x,datlen);plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y

30、)title(过滤后语音信号频域波形);4.4.4 信号经过高通 fir%窗函数法设计窗函数法设计高通滤波器（高通fir）clear all;x1,fs,bits=wavread(c:matlab7workcricket.wav);datlen=length(x1);sound(x1,fs,bits);y1=fftshift(fft(x1);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y1);title(原始语音信号频谱)n=le

31、ngth(x1)-1;t=0:1/fs:n/fs;fc=4800;fp=5000;as=100;ap=1;fs=22050; %采样频率必须变大（n-1）ws=2*pi*fc/fs; wp=2*pi*fp/fs;bt=wp-ws;n=ceil(8*pi/bt);wc=(wp+ws)/2/pi;b,a=fir1(n,wc,high);figure(2);freqz(b,a);x=fftfilt(b,x1);sound(x,fs,bits);figure(3) %画出滤波后的信号subplot(2,1,1)plot(x);title(过滤后语音信号时域波形);subplot(2,1,2)y=fftshift(fft(x,datlen);plot(-datlen/2:datlen/2-1)/datlen*fs,abs(y);title(过滤后语音信号频域波形);4.5 撰写5000-8000字课程设计报告5. 总结 首先，对一段语音信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。 通过数