基于MATLAB语音信号采集课程设计

1、数字信号处理课程设计报告课设题目：语音信号的采集与处理 学 院：信息科学与工程学院 专 业：电子信息工程 班 级：电子0903 姓 名： 学 号： 指导教师：张海云、李锦彬 时 间：2012.12.312013.01.11一.引言介绍.二. matlab工具简介.三. 课程设计任务.3.1语音信号的采集.3.2语音信号的频谱分析.3.3设计数字滤波器和画出其频率响应给出各滤波器的性能指标.3.4用滤波器对信号进行滤波.3.5回放语音信号，分析滤波前后的语音变化.四. 课程设计原理及方案.4.1用窗函数法设计fir滤波器.4.2用双线性变换法设计iir数字滤波器.五. 语音信号的采集及频谱分析.

2、六. 设计滤波器及对语音信号进行滤波.6.1 fir低通滤波器设计及信号滤波.6.2fir高通滤波器设计及信号滤波.6.3 fir带通滤波器设计及信号滤波.6.4 iir低通滤波器设计及信号滤波.6.5 iir高通滤波器设计及信号滤波.6.6 iir带通滤波器设计及信号滤波.七. 课程设计心得体会.八. 参考文献.一引言随着信息时代和数字世界的到来，数字信号处理已成为今一门极其重要的学科和技术领域。数字信号处理在通信语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理应用中，数字滤波器十分重要并已获得广泛应用。二matlab工具简介1matlab

3、是矩阵实验室（matrixlaboratory）之意，现已发展成为适合多学科，多种工作平台的功能强大的大型软件，已经成为线性代数，自动控制理论，数理统计，数字信号处理，时间序列分析，动态系统仿真等高级课程的基本教学工具；2matlab的语言特点（1）语言简洁紧凑，使用方便灵活，库函数极其丰富。matlab程序书写形式自由，利用起丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。（2）运算符丰富。由于matlab是用c语言编写的，matlab提供了和c语言几乎一样多的运算符，灵活使用matlab的运算符将使程序变得极为简短

4、。 （3）matlab既具有结构化的控制语句（如for循环，while循环，break语句和if语句），又有面向对象编程的特性。 （4）程序限制不严格，程序设计自由度大。例如，在matlab里，用户无需对矩阵预定义就可使用。 （5）程序的可移植性很好，基本上不做修改就可在各种型号的计算机和操作系统上运行。 （6）matlab的图形功能强大。在fortran和c语言里，绘图都很不容易，但在matlab里，数据的可视化非常简单。matlab还具有较强的编辑图形界面的能力。 （7）matlab的缺点是，它和其他高级程序相比，程序的执行速度较慢。由于matlab的程序不用编译等预处理，也不生成可执行文

5、件，程序为解释执行，所以速度较慢。 （8）功能强大的工具箱是matlab的另一特色。matlab包含两个部分：核心部分和各种可选的工具箱。核心部分中有数百个核心内部函数。其工具箱又分为两类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能，图示建模仿真功能，文字处理功能以及与硬件实时交互功能。功能性工具箱用于多种学科。而学科性工具箱是专业性比较强的。（9）源程序的开放性。开放性也许是matlab最受人们欢迎的特点。除内部函数以外，所有matlab的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。三课程设计任务3.1语音信号的

6、采集 利用windows下的录音机，录制一段自己的话音，时间在3秒左右，然后在matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数，采样率8000hz，单通道，格式*.wav。注：如果没有办法录制自己语音，也可以用现成的语音也可以，文件格式 *.wav3.2语音信号的频谱分析 在matlab中，可以利用函数fft对信号进行快速傅立叶变换，得到信号的频谱特性，要求学生首先画出语音信号的时域波形，然后对语音信号进行频谱分析。3.3设计数字滤波器和画出其频率响应给出各滤波器的性能指标； 给定滤波器的性能指标如下： (1)低通滤波器的性能指标：fb=1000hz,f

7、c=1200hz,as=100db,ap=1db. (2)高通滤波器的性能指标：fb=2700hz, fc=3000hz,as=100db,ap=1db. (3)带通滤波器的性能指标： fb1=1200hz, fb2=3000hz, as=100db,ap=1db. fc1=1000hz, fc2=3200hz,as=100db, ap=1db. 采用窗函数法和双线性变换法设计上面要求的3种滤波器，并画出滤波器的频率响应； 3.4用滤波器对信号进行滤波 然后用自己设计的滤波器对采集到的信号进行滤波，画出滤波后信号的时域波形及频谱，并对滤波前后的信号进行对比，分析信号的变化； 3.5回放语音信号

8、，分析滤波前后的语音变化； 四滤波器设计原理及设计方案4.1用窗函数法设计fir滤波器根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度n（或阶数m=n-1），窗函数类型可根据最小阻带衰减as独立选择，因为窗口长度n对最小阻带衰减as没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度n，设待求滤波器的过渡带宽为w，它与窗口长度n近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的n，在n和窗函数类型确定后，即可调用matlab中的窗函数求出窗函

9、数wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2采样n点，采用离散傅里叶反变换（idft）即可求出。用窗函数wd(n)将hd(n)截断，并进行加权处理，得到如果要求线性相位特性， 则h(n)还必须满足：根据上式中的正、 负号和长度n的奇偶性又将线性相位fir滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(n-1-n)一类，而不能选h(n)=

10、-h(n-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。4.2用双线性变换法设计iir数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从s平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/t/t之间，再用z=est转换到z平面上。也就是说，第一步先将整个s平面压缩映射到s1平面的-/t/t一条横带里；第二步再通过标准变换关系z=es1t将此横带变换到整个z平面上去。这样就使s平面与z平

11、面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1所示。图1双线性变换的映射关系为了将s平面的整个虚轴j压缩到s1平面j1轴上的-/t到/t段上，可以通过以下的正切变换实现 （1）式中,t仍是采样间隔。当1由-/t经过0变化到/t时，由-经过0变化到+，也即映射了整个j轴。将式（1）写成将此关系解析延拓到整个s平面和s1平面，令j=s，j1=s1，则得再将s1平面通过以下标准变换关系映射到z平面z=es1t从而得到s平面和z平面的单值映射关系为： (2) (3)式（2）与式（3）是s平面与z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变

12、换式（1）与式（2）的双线性变换符合映射变换应满足的两点要求。首先,把z=ej，可得 (4)即s平面的虚轴映射到z平面的单位圆。其次，将s=+j代入式（4），得因此由此看出，当0时，|z|0时，|z|1。也就是说，s平面的左半平面映射到z平面的单位圆内，s平面的右半平面映射到z平面的单位圆外，s平面的虚轴映射到z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点:双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为s平面与z平面是单值的一一对应关系。s平面整个j轴单值地对应于z平面单位圆一周，即频率轴是单值变换关

13、系。这个关系如式（4）所示，重写如下：上式表明，s平面上与z平面的成非线性的正切关系，如图2所示。由图2看出，在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但当进一步增加时，增长得越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。图2双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（4）及图2所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次

14、，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3所示。图3双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。五语音信号的采集及频谱分析用windows下的录音机，用单声道录制一段音乐或声音，

15、时间在5s内。然后matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。对语音信号进行快速傅立叶变换，在一个窗口同时画出信号的时域波形图和频谱图，分析语音信号的频谱特点y,fs,nbits=wavread(music.wav);%加载语音信号利用函数wavread对语音信号进行采样wavplay(y,fs);%回放语音信号n=length(y);%语音信号长度y=fft(y,n);%用傅里叶变换求频域subplot(2,1,1);plot(y);title(图1 原始信号时域波形图);subplot(2,1,2);plot(abs(y);title(图2 原

16、始信号频谱图);六.设计滤波器及对语音信号进行滤波6.1fir低通滤波器设计及信号滤波x1,fs,bits=wavread(music.wav); derta_fs = fs/length(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致fs=fs; fp1=1000; fs1=1200; as1=100;wp1=2*pi*fp1/fs; %ws1=2*pi*fs1/fs; %bf1=ws1-wp1;wc1=(wp1+ws1)/2;m1=ceil(as1-7.95)/(2.286*bf1)+1;%按凯泽窗计算滤波器阶数n1=m1+1;beta1=0.1102*(as1

17、-8.7); window=(kaiser(n1,beta1); %求凯泽窗窗函数b1=fir1(m1,wc1/pi,window);% wc1/pi为归一化,窗函数法设计函数figure(2); freqz(b1,1,512); %h,w=freqz(b,a,n),(1)中b和a分别为离散系统的系统函数分子、分母多项式的系数向量，返回量h则包含了离散系统频响在 0pi范围内n个频率等分点的值（其中n为正整数），w则包含了范围内n个频率等分点。调用默认的n时，其值是512。title(fir低通滤波器的频率响应); x1_low = filter(b1,1, x1);%对信号进行低通滤波 ,y

18、 = filter(b,a,x) ，输入x为滤波前序列，y为滤波结果序列，b/a 提供滤波器系数，b为分子， a为分母 sound(x1_low,fs,bits); figure(3);subplot(2,1,1);plot(x1_low);title(信号经过fir低通滤波器(时域); subplot(2,1,2);plot(-fs/2:derta_fs: fs/2-derta_fs,abs(fftshift(fft(x1_low); title(信号经过fir低通滤波器（频域）);6.2 fir高通滤波器设计及信号滤波x1,fs,bits=wavread(music.wav); derta

19、_fs = fs/length(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致fs=fs; as2=100;fp2=3000; fs2=2700;wp2=2*pi*fp2/fs; ws2=2*pi*fs2/fs; bf2=wp2-ws2;wc2=(wp2+ws2)/2;m2=ceil(as2-7.95)/(2.286*bf2)+1;%按凯泽窗计算滤波器阶数n2=m2+1;beta2=0.1102*(as2-8.7); window=(kaiser(n2,beta2); %求凯泽窗窗函数b2=fir1(m2,wc2/pi,high,window);figure(4);

20、 freqz(b2,1,512);%数字滤波器频率响应 title(fir高通滤波器的频率响应); x1_high = filter(b2,1,x1);%对信号进行高通滤波sound(x1_high,fs,bits);figure(5);subplot(211);plot(x1_high);title(信号经过fir高通滤波器(时域);subplot(212);plot(-fs/2:derta_fs: fs/2-derta_fs,abs(fftshift(fft(x1_high);title(信号经过fir高通滤波器（频域）);6.3 fir带通滤波器设计及信号滤波fs=fs;as3=100;

21、fp3=1200,3000;fs3=1000,3200;wp3=2*pi*fp3/fs; ws3=2*pi*fs3/fs; bf3=wp3(1)-ws3(1);wc3=wp3+bf3/2;m3=ceil(as3-7.95)/(2.286*bf3)+1;%按凯泽窗计算滤波器阶数n3=m3+1;beta3=0.1102*(as3-8.7); window=(kaiser(n3,beta3); %求凯泽窗窗函数b3=fir1(m3,wc3/pi,bandpass,window);%带通滤波器 figure(6);freqz(b3,1,512);%数字滤波器频率响应 title(fir带通滤波器的频率

22、响应);x1_daitong = filter(b3,1,x1);%对信号进行带通滤波sound(x1_daitong,fs,bits);figure(7);subplot(211);plot(x1_daitong);title(信号经过fir带通滤波器(时域); subplot(212);plot(-fs/2:derta_fs: fs/2-derta_fs,abs(fftshift(fft(x1_daitong);title(信号经过fir带通滤波器（频域）);6.4 iir低通滤波器设计及信号滤波x1,fs,bits=wavread(music.wav); derta_fs = fs/le

23、ngth(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致fs=fs;fp1i=1000; fs1i=1200;wp1i=2*pi*fp1i/fs; ws1i=2*pi*fs1i/fs; rp1i=1; rs1i=100; ts=1/fs; wp1i=2/ts*tan(wp1i/2);ws1i=2/ts*tan(ws1i/2); %按频率转换公式进行转换，预畸变 n1i,wn1i=cheb1ord(wp1i,ws1i,rp1i,rs1i,s); %计算模拟滤波器的最小阶数 b1i,a1i=cheby1(n1i,rp1i,wn1i,s);%设计模拟原型滤波器 bz1i,

24、az1i=bilinear(b1i,a1i,fs); %运用双线性变换法得到数字滤波器传递函数 figure(10);freqz(bz1i,az1i,512,fs); title(切比雪夫1型低通滤波器的频率响应);x1_lowi = filter(bz1i,az1i, x1);%对信号进行低通滤波 ,y = filter(b,a,x) ，输入x为滤波前序列，y为滤波结果序列，b/a 提供滤波器系数，b为分子， a为分母 sound(x1_lowi,fs,bits); figure(11); subplot(211);plot(x1_lowi);title(信号经过iir低通滤波器(时域);

25、subplot(212);plot(-fs/2:derta_fs: fs/2-derta_fs,abs(fftshift(fft(x1_lowi); title(信号经过iir低通滤波器（频域）);6.5 iir高通滤波器设计及信号滤波待添加的隐藏文字内容1x1,fs,bits=wavread(music.wav); derta_fs = fs/length(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致% -iir高通按钮消息响应函数.fs=fs;fp2i=3000;%通带截止频率fs2i=2700;%组带截止频率rp2i=1; rs2i=100;wp2i2=2*

26、pi*fp2i/fs;ws2i2=2*pi*fs2i/fs;wp2i=tan(wp2i2/2);ws2i=tan(ws2i2/2);wphi=1 ; %归一化wshi=wp2i/ws2i;%高通频率转换成低通频率指标%估计滤波器的阶数n2i,wn2i=cheb1ord(wphi,wshi,rp2i,rs2i,s);%设计滤波器b2i,a2i=cheby1(n2i,rp2i,wn2i,s);numi2,deni2=lp2hp(b2i,a2i,wp2i);bz2i,az2i=bilinear(numi2,deni2,0.5)figure(12);freqz(bz2i,az2i,512,fs); t

27、itle(切比雪夫1型iir高通滤波器的频率响应);x1_highi = filter(bz2i,az2i, x1);%对信号进行高通滤波 ,y = filter(b,a,x) ，输入x为滤波前序列，y为滤波结果序列，b/a 提供滤波器系数，b为分子， a为分母 sound(x1_highi,fs,bits); figure(13); subplot(211);plot(x1_highi);title(信号经过iir高通滤波器(时域); subplot(212);plot(-fs/2:derta_fs: fs/2-derta_fs,abs(fftshift(fft(x1_highi); tit

28、le(信号经过iir高通滤波器（频域）);6.6 iir带通滤波器设计及信号滤波x1,fs,bits=wavread(music.wav); derta_fs = fs/length(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致fs=fs;fdp1i=1200;fdp2i=3000;%通带截止频率fds1i=1000;fds2i=3200;%组带截止频率rp3i=1; rs3i=100; wdp1i=2*pi*fdp1i/fs;wdp2i=2*pi*fdp2i/fs; wds2i=2*pi*fds1i/fs;wds2i=2*pi*fds2i/fs;wdp1i=ta

29、n(wdp1i/2);wdp2i=tan(wdp2i/2);wds1i=tan(wds2i/2); wds2i=tan(wds2i/2); %按频率转换公式进行转换，预畸变 b=wdp2i-wdp1i; w3i = wdp1i*wdp2i/wds1i; %调整截止频率参数。 w0 = (wdp1i*wdp2i)0.5; wp3i = 1; %归一化，求低通原型的参数。ws3i = -(wdp1i*wdp2i)-wds2i2)/(b*wds2i); n3i,wn3i=buttord(wp3i,ws3i,rp3i,rs3i,s); %计算模拟滤波器的最小阶数 b3i,a3i=butter(n3i,

30、wn3i,s);%设计模拟原型滤波器 numi3,deni3= lp2bp(b3i,a3i,w0,b); bz3i,az3i = bilinear(numi3,deni3,0.5); %双线性变换 h,w = freqz(bz3i,az3i,512); figure(14);subplot(211); plot(w/pi,20*log10(abs(h);grid; ylabel(gain); title(巴特沃兹iir带通滤波器) axis(0 1 -80 5); subplot(212); plot(w/pi,angle(h);grid; xlabel(w/pi),ylabel(phase)

31、;x1_daitongi = filter(bz3i,az3i, x1);%对信号进行带通滤波 ,y = filter(b,a,x) ，输入x为滤波前序列，y为滤波结果序列，b/a 提供滤波器系数，b为分子， a为分母 sound(x1_daitongi,fs,bits); figure(15); subplot(211); plot(x1_daitongi);title(信号经过iir带通滤波器(时域); subplot(212); plot(-fs/2:derta_fs: fs/2-derta_fs,abs(fftshift(fft(x1_daitongi); title(信号经过iir带通滤波器（频域）)