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

1、长沙学院电子与通信工程系- -题 目系（部）专业（班级）姓 名学 号指导老师起止日期数字信号处理课程设计说明书数字信号处理课程设计电子与通信工程系10级通信2班周斌2010043205李广柱、刘光灿、陈威兵黄飞江、张刚林、冯璐2013.4.15-2013426目录 TOC o 1-5 h z 绪论 -4-设计作品名称 -4-课程设计内容 -4-3.1设计思想 -4- HYPERLINK l bookmark6 o Current Document 3.2设计原理 -5- HYPERLINK l bookmark8 o Current Document 3.2.1无限脉冲响应数字滤波器设计 -5

2、-3.2.1.1巴特沃斯滤波器设计： -5-巴特沃斯滤波器的特性： -5-双线性变换法： -7- HYPERLINK l bookmark45 o Current Document 3.2.2有限脉冲响应数字滤波器设计 -9- HYPERLINK l bookmark65 o Current Document 课程设计设计步骤及结果分析 -12- HYPERLINK l bookmark67 o Current Document 4.1语音信号的采集及其时域、频域分析 -12-程序 -12 -图像 -12- HYPERLINK l bookmark83 o Current Document 4

3、.2设计数字滤波器和画出其频率相映 -13- HYPERLINK l bookmark85 o Current Document 无限低通数字滤波器 -13 -程序 -13 -图像 -14 - HYPERLINK l bookmark89 o Current Document 无限高通数字滤波器 -14 -程序 -14 -图像 -15 - HYPERLINK l bookmark93 o Current Document 有限低通数字滤波器 -16 -程序 -16 -图像 -16 - HYPERLINK l bookmark97 o Current Document 有限高通数字滤波器 -17

4、-程序 -17 -图像 -17 - HYPERLINK l bookmark99 o Current Document 4.3用滤波器对信号进行滤波并比较前后波形 -18 -经过低通 IIR - 18 -经过高通 IIR - 19 -经过低通 FIR - 20 -经过高通 FIR - 21 -编制实现上述任务的相应的总程序 -21 -信号经过低通IIR - 21 -信号经过高通IIR - 22 -信号经过低通FIR - 24 -信号经过高通 FIR - 25 - TOC o 1-5 h z HYPERLINK l bookmark117 o Current Document 4.5撰写5000

5、-8000字课程设计报告 -26 -总结 -26 - HYPERLINK l bookmark121 o Current Document 存在的建议及不足 -27 - HYPERLINK l bookmark123 o Current Document 参考文献 -27 -摘要对一段语音信号进行采样；画出采样后语音信号的时域波形和频谱图； 给定 滤波器的性能指标，采用窗函数法和双线性变换设计滤波器，并画出滤波器的频 率响应；然后用自己设计的滤波器对采集的信号进行滤波， 画出滤波后信号的时 域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号； 最后，设计一个信号处理系统界面

6、。通过数字信号处理的课程设计，使学生对信号的采集，处理，传输，显示， 存储和分析等有一个系统的掌握和理解。巩固和运用数字信号处理课程中的理论 知识和实验技能，掌握最基本的数字信号处理的理论和方法， 培养学生发现问题， 分析问题和解决问题的能力。为以后的工作和学习打下基础。本课题主要基于 MATLA平台，实现信号分析与滤波器的设计。并改变滤波器参数或特性 （低通、 高通），实现不同的滤波要求。自己录播一段语音信号，对其进行采样，进行频 谱分析，设计有限和无限两种不同的数字滤波器对信号进行滤波处理，分析所设计滤波器（画出了频率特性曲线），并对信号进行滤波，画出滤波后信号的时域波 形和频谱，并对滤波

7、前后的信号进行对比，分析信号的变化，分析的结果表明初 步实现了设计目标。关键词：无限脉冲响应数字滤波器、有限脉冲响应数字滤波器、巴特沃斯滤波 器、窗函数法、双线性法1緒论对一段语音信号进行采样； 画出采样后语音信号的时域波形和频谱图； 给定滤波器的性 能指标，采用窗函数法和双线性变换设计滤波器， 并画出滤波器的频率响应； 然后用自己设 计的滤波器对采集的信号进行滤波， 画出滤波后信号的时域波形和频谱， 并对滤波前后的信 号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。通过数字信号处理的课程设计，使学生对信号的采集，处理，传输，显示，存储和分析等有一个系统的掌握和理解。

8、巩固和运用数字信号处理课程中的理论知识和实验技能，掌握最基本的数字信号处理的理论和方法，培养学生发现问题， 分析问题和解决问题的能力。为以后的工作和学习打下基础。本课题主要基于MATLAB台，实现信号分析与滤波器的设计。 并改变滤波器参数或特性（低通、高通），实现不同的滤波要求。 自己录播一段语音信号，对其进行采样，进行频谱分析，设计有限和无限两种不同的数字滤波器对信号进行滤波处理， 分析所设计滤波器（画出了频率特性曲线），并对信号进行滤波，画出滤波后信号的时域波形 和频谱，并对滤波前后的信号进行对比，分析信号的变化，分析的结果表明初步实现了设计目标。2设计作品名称语音信号的采集、分析与处理3

9、课程设计内容3.1设计思想采样一段语音信号； 画出语音信号的时域波形和频谱图； 给定滤波器的性能指标， 设计 数字滤波器，并画出滤波器的频率响应； 然后用设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱。关键技术：频谱图的理解；设计数字滤波器；数字滤波的方法；解决思路：对语音号进行快速傅里叶变换，得到信号的频谱特性；在MATLAB境中可以利用函数firl设计FIR滤波器，可以利用函数butter，chebyl和ellip 设计IIR滤波器；利 用MATLAB中的函数freqz画出各滤波器的频率响应。3.2设计原理321无限脉冲响应数字滤波器设计3.2.1.1巴特沃斯滤波器设计:

10、巴特沃斯滤波器的特点是同频带的频率响应去想最平滑，如下图所示图1 一级巴特沃斯低通滤波器的波得图2ikl dcr .WcKkr 4tholcr Sih enter0(11图2 一级至五级巴特沃斯低通滤波器3.1.1.2巴特沃斯滤波器的特性：没有起伏，而在组频带则巴特沃斯滤波器的特点是同频带内的频率响应曲线最为平坦,逐渐下降为零。在振幅的对数对角频率的波特图上，从某一边界见频率开始，振幅随着角频率的增加而逐渐减少，趋向于负无穷大。一阶巴特沃斯滤波器的衰减率为每倍频20分贝，二阶巴特沃斯滤波器的衰减率为每倍频12分贝，三阶的衰减率为每分贝 18分贝，如此类推，巴特沃斯滤波器的振幅对角频率单 调下降

11、，并且滤波器的结束越高，在组频带振幅衰减速度越快，其他滤波器高阶的振幅对角 频率图和低阶数的振幅对角频率有不同的形状。上述函数的特点是等距离分布在半径为I 的圆上。因此，极点用下式表示为匚 j(2k 1)Skk =0,1,2, N -1Ha(s)的表示式:onHa(s)f丨丨(s - sk )k=0为了使设计公式和图表统一，将频率归一化。巴特沃斯滤波器采用3dB截止频率xc归一化，归一化后的系统函数为Ga (j ) = n 41 sc (S J k4c A令p = j = s.黒， -y, 称为归一化频率，p称为归一化复变量，这样巴特沃斯滤波器的归一化低通原型系统函数为1Ga N -4(p -

12、 Pk)k=0式中，Pk二s ic，为归一化极点，用下式表示j 7(2 Pk =e2k 1)2Nk =0,1,2, N -1然后将巴特沃斯滤波器转换成高通滤波器，转换公式如下:P =，p： ； ph / S在p平面与s平面虚轴上的频率关系为：其中，11 ph为希望设计的高通滤波器H HP(S)的通带边界频率。H hp (s) = G( p) pGm/s3.1.1.3双线性变换法:双线性变换法是将 S平面压缩变换到某一中介S1平面的一条横带里，再通过标准变换关系z二exp(s1*T)将此带变换到整个 z平面上去，这样就使s平面与z平面之间建立对应的单值关系，消除了多值变换性。为了将s平面的j-1

13、轴压缩到s1平面的j1轴上的一 pi T到pi T 一段上，可以通过以下的正切变换来实现： TOC o 1-5 h z HYPERLINK l bookmark49 o Current Document 21tan ( JT)T2这样当门1由- piT经0变化到pi T时，门由-经过0变化到，也映射到了整个s1 T2 1 -es1 TT 1 ej轴。将这个关系延拓到整个s平面和s1平面，则可以得到2T HYPERLINK l bookmark31 o Current Document stan (s1)=T2再将s1平面通过标准变换关系映射到z平面，即令z =exp(s1* T)得到同样对z求

14、解，得到2 sZ2 s T这样的变换叫做双线性变换。为了验证这种映射具有 s平面的虚轴映射到 Z平面单位圆上的特性，考虑 s = ji】，Z ，得j11门=-ta n丄T 2除了使s平面的虚轴映射到单位圆上之外,s平面的左半部分映射到单位圆的内部,s平面的右半部分映射到单位圆的外部。如图所示图3双线性变化映射关系示意图观察式子才，发现s的实部为负时，因子的幅度小于1,相当于单位圆的内部。反之，当s的实部为负时，该比值的幅度大于1，相当于单位圆的外部。这样就可以看出使用双线性变换可从稳定的模拟滤波器得到稳定的数字滤波器。双线性变换法还避免了使用脉冲响应不变法所遇到的混叠问题，因为它把s平面的这个

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

16、要使s平面与z平面建立 对 应的单值关系，即求出 s = f(z)，然后将其代入 G(s)就可以求得H(z)，即H(z)=G(s)s=g3.2.2有限脉冲响应数字滤波器设计如果所希望的滤波器的理想的频率响应函数为Hd e，则其对应的单位脉冲响应为hd n H d ej ej d (4.1)2窗函数设计法的基本原理是用有限长单位脉冲响应序列h n逼近hd n。由于hd n往往是无限长序列，而且是非因果的，所以用窗函数n将hd n截断，并进行加权处理，得到：(4.2)h n = hd n nhn就作为实际设计的 FIR数字滤波器的单位脉冲响应序列，其频率响应函数 Hej为N _1H e，八 h n

17、ej n(4.3 )n=0式中，N为所选窗函数n的长度。我们知道，用窗函数法设计的滤波器性能取决于窗函数n的类型及窗口长度 N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。表1各种窗函数的基本参数窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134 n /N-12三角形窗-258n /N-25汉宁窗-318n /N-44哈明窗-418n /N-53不莱克曼窗-5712n /N-74凯塞窗(a =7.865)-5710n /N-80这样选定窗函数类型和长度 N之后，求出单位脉冲响

18、应 h n二hd n n ,并按照式 (4.3 )求出He。 H e是否满足要求，要进行演算。一般在h n尾部加零使长度满足2的整数次幕，以便用FFT计算H e。如果要观察细节，补零点数增多即可。如果H ej 不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。如果要求线性相位特性，则 h n还必须满足h n = h N -1 - n根据上式中的正、负号和长度N的奇偶性又将线性相位 FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择h n i；=h N -1 - n这一类，而不能选择 h n二_h N -1 - n这一类。主程序

19、框图如图4.1所示。其中幅度特性要求用dB表示。图5主程序框图H (k) =DFTh(n)H(k)二HR(k) jHi(k)H(k) =#HR(k) +H：(k)一i一、2-画图时，用20lgH(k)打印幅度特性。第 k点对应的频率 纵=k。为使曲线包络更接N近He的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h n末尾补零的方法，使长度变为64,再进行64点DFT，则可以得到更精确的幅度衰减特性曲线。4.课程设计设计步骤及结果分析4.1语音信号的采集及其时域、频域分析对一段语音信号，时间在 2s内。然后在MATLAB软件平台下，利用函数wavread对语 音信号进行

20、采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解掌握采样频率、采样位数等概念。要求首先画出采样后语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。4.1.1程序x1,fs,bits=wavread(C:MATLAB7workcricket.wav);% 原信号sou nd(x1,fs,bits); %调用sou nd函数播放原始语音figure(1);plot(x1);%做原始语音信号的时域图形title(原始语音信号时域图形)；xlabel(时间 t);ylabel(音量 n);figure(2);y仁fft(x1);% 做

21、 length(x1)点的 FFTy仁fftshift(yl); %平務 卫频邛门心人oderta_fs = fs/length(x1); %设置频谱的间隔 分辨率plot(-fs/2:derta_fs: fs/2-derta_fs,abs(y1);% 画出原始语音信号的频谱图title(原始语音信号的频谱图)；grid on;4.1.2图像原始语音信言时域波形600T4J0200-1原始语音信号频谱|0.511.5x 10*4.2设计数字滤波器和画出其频率相映给出各滤波器的性能指标： 低通滤波器性能指标fb = 1000Hz, fc = 1200Hz, As = 100dB, Ap = 1d

22、B。(2)高通滤波器性能指标fc = 4800Hz , fb = 5000Hz, As = 100dB, Ap = 1dB。要求用窗函数法和双线性变换法设计上面要求的2种滤波器。在MATLAB中，可以利用函数fir1设计FIR滤波器，可以利用函数 butter设计IIR滤波器；利用 MATLAB中的函 数freqz画出各滤波器的频率响应。4.2.1无限低通数字滤波器421.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,k);bl,

23、al=lp2lp(bap,aap,wc);b,a=bili near(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(IIR低通滤波器频率响应);421.2图像IIR低通滤波器频率晌应无限高通数字滤波器程序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=bili

24、 near(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(IIR高通滤波器频率响应);4.222图像职咼通滤液器频率响应141 210.8060.40.20020004000600080001000012000423有限低通数字滤波器423.1程序N=le ngth(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.

25、7); %alph用来控制窗的形状N=ceil(As-8)/2.285/Bt); %凯塞窗滤波器阶数 N,(h(n)的长度为N+1)wn= kaiser(N+1,alph);ws=(wp+wc)/2/pi;%计算通带截止频率(关于 pi归一化)hn=fir1(N,ws,w n);figure(2)%画出滤波器频谱图freqz(h n,1)423.2图像111111111111|1111FliiiIIi|i1144V11114!1111li1111屮11i|11*1：LIFII1)1)|iPh|i|iVji|i11l!bw-訂I111Vi 1!111-1!|1111i1iI11-i1i1ii1i

26、iii!iLjiiiIl11I1i|11i1r T *|V|i1!hh|i|iVVfliI114F |-I1111-|i|ii|ijii0.10.20.30.40.50.60.70.30.91001Normalized Frequency (x rad sample)OJ 020.3040 50 60.70 80 91NormaHzed Frequency (x rad/sample)o o o JI2030-424有限高通数字滤波器424.1程序N=le ngth(x1)-1;t=O:1/fs:N/fs;fc=4800;fp=5000;As=100;Ap=1;Fs=22050; %采样频率必

27、须变大(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);424.2图像0 20.30 40 50 6070 80 91iNormaiized Frequency (x?r rad/sample(忻 x-sflJLICI0.20.30.40 50.60.7 Q 80 91Normalized Frequency rad/ssmple4.3用滤波器对信号进行滤波并比较前后波形要求用自己设计的各种滤波器分别对采集的信号

28、进行滤波，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。经过低通IIR过滤后语音信号时域波形过滤后语音信号频域波彫经过高通IIR过滤石语音信号时域波形00 511.622.533.544 5x 10彳 5J-0.500.511.5 x 10过滤后语音信号频域波形Ji4.3.3 经过低通FIR过滤后语音信号时域波形0 50-0 5-111|L|L|I|lII1llIJIillI00 511,522 533 544x 104 过滤后语音信号频域波形1 1|11 1600400200500.5x1044.3.4 经过高通FIR过滤

29、后语音信号时域淤形过滤后语音信号频域渡形x 1044.4编制实现上述任务的相应的总程序4.4.1信号经过低通IIR9无限脉冲响应数字低通滤波器(低通 IIR)clc;clear all;x1,fs,bits=wavread(C:MATLAB7workcricket.wav);datle n=len gth(x1);soun d(x1,fs,bits);y1=fftshift(fft(x1);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datle n/2:datle n/2-1)/datle n*fs

30、,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);bap,aap=zp2tf(z,p,k);bl,al=lp2lp(bap,aap,wc);b,a=bili near(bl,al,fs);h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(IIR低通滤波器频率响应);X=filter(b,a,x1);%将信号经过滤波器sou nd(X,fs,bits);figure(3)%画出

31、滤波后的信号subplot(2,1,1)plot(X);title(过滤后语音信号时域波形);subplot(2,1,2)y3=fftshift(fft(X)plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y3);title(过滤后语音信号频域波形);4.4.2信号经过高通IIR溉限脉冲响应数字高通滤波器(高通IIR)clear all;x1,fs,bits=wavread(C:MATLAB7workcricket.wav);datle n=len gth(x1);soun d(x1,fs,bits);y仁fftshift(fft(x1);figure(

32、1)subplot(2,1,1) plot(xl);title(原始语音信号时域波形);subplot(2,1,2)plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y1); title(原始语音信号频谱)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=bili near(bl,al,fs);

33、h,w=freqz(b,a);figure(2)plot(w*fs/(2*pi),abs(h)title(IIR 高通滤波器频率响应);X=filter(b,a,x1);%将信号经过滤波器sou nd(X,fs,bits);figure(3)%画出滤波后的信号subplot(2,1,1) plot(X);title（过滤后语音信号时域波形）;subplot(2,1,2)y3=fftshift(fft(X) plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y3);title(过滤后语音信号频域波形);4.4.3信号经过低通FIR%言号经过用窗函数法设计低通

34、FIRclc;clear all;x1,fs,bits=wavread(C:MATLAB7workcricket.wav);datle n=len gth(x1);soun d(x1,fs,bits);y1=fftshift(fft(x1);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y1);title(原始语音信号频谱)N=le ngth(x1)-1;t=0:1/fs:N/fs;fp=1000;fc=1500;As=100

35、;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); % wn= kaiser(N+1,alph);ws=(wp+wc)/2/pi; %hn=fir1(N,ws,w n);figure(2)%freqz(h n,1)X=co nv(hn ,x1);%sou nd(X,fs,bits);figure(3)%subplot(2,1,1) plot(X);计算过渡带宽用来控制窗的形状凯塞窗滤波器阶数 N,(h(n)的长度为N+1)计

36、算通带截止频率(关于pi归一化)画出滤波器频谱图将信号经过滤波器画出滤波后的信号title(过滤后语音信号时域波形);subplot(2,1,2)y=fftshift(fft(X,datle n);plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y)title(过滤后语音信号频域波形);4.4.4信号经过高通FIR9窗函数法设计窗函数法设计高通滤波器(高通 FIR)clear all;x1,fs,bits=wavread(C:MATLAB7workcricket.wav);datle n=len gth(x1);soun d(x1,fs,bits);y1

37、=fftshift(fft(x1);figure(1)subplot(2,1,1)plot(x1);title(原始语音信号时域波形);subplot(2,1,2)plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y1);title(原始语音信号频谱)N=le ngth(x1)-1;t=0:1/fs:N/fs;N-1)fc=4800;fp=5000;As=100;Ap=1;Fs=22050; %采样频率必须变大(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);soun d(x,fs,bits);figure(3)%画出滤波后的信号subplot(2,1,1)plot(x);title(过滤后语音信号时域波形);subplot(2,1,2)y=fftshift(fft(x,datle n);plot(-datle n/2:datle n/2-1)/datle n*fs,abs(y);