数字信号处理课程设计语音信号的频谱分析

《数字信号处理课程设计语音信号的频谱分析》由会员分享，可在线阅读，更多相关《数字信号处理课程设计语音信号的频谱分析（21页珍藏版）》请在装配图网上搜索。

1、目 录1.课程设计的目的12.课程设计的要求1（1）熟悉离散信号和系统的时域特性1（2）掌握序列快速傅里叶变换FFT方法1（3）学会MATLAB的使用，掌握MATLAB的程序设计方法1（4）利用MATLAB对语音信号进行频谱分析1（5）掌握MATLAB设计FIR和IIR数字滤波器的方法13.课程设计报告内容14.设计步骤1（1）语音信号的获取1（2）语音信号的频谱分析3（3）设计FIR和IIR数字滤波器并画出频率响应4（4）对信号进行滤波8（5）被污染语音信号的频谱分析13（6）回放语音信号135.调试分析135.1等间隔采样程序的编写135.2 .m文件的命名145.3标点符号的使用155.

2、4滤波器参数的使用156.结果分析与体会166.1总结166.2不足176.3思考176.4结束语17参考文献18附录：1920语音信号的频谱分析1.课程设计的目的综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解。2.课程设计的要求（1）熟悉离散信号和系统的时域特性。（2）掌握序列快速傅里叶变换FFT方法。（3）学会MATLAB的使用，掌握MATLAB的程序设计方法。（4）利用MATLAB对语音信号进行频谱分析。（5）掌握MATLAB设计FIR和IIR数字滤波器的方法。3.课程设计报告内容录

3、制一段自己的语音信号，取不同的数据点对语音信号进行频谱分析；对所有数据进行插值和抽取处理，改变抽样率再对信号进行频谱分析；设计FIR和IIR数字滤波器，并对被抽样后的语音信号进行滤波，分析滤波后信号的时域和频域特征，回放语音信号。4.设计步骤（1）语音信号的获取用Windows下录音机录制自己的声音，发音为汉语韵母a，保存在D盘命名为5.wav。通过使用wavread函数读取语音信号，实现语句如下：x,fs,bits=wavread(d:5.wav);为了得到原始信号的波形和频谱，可运行如下语句：xx=fft(x);subplot(2,1,1);plot(x);title(原始信号波形);su

4、bplot(2,1,2);plot(abs(xx); title(原始信号频谱);运行后得到的波形、频谱图如下：原始的语音信号经过wavread函数的处理，由连续信号变成离散信号，抽样率如图中fs所示为22050 Hz/s，为了使抽样率满足8000 Hz/s的要求，需要用到resample函数，实现语句如下：y=resample(x,8000,22050);经过重采样，原信号的抽样率变成8000Hz/s。由于采样频率的改变，在模拟信号中选取的点数也有所不同,采样频率由22050 Hz/s变成8000 Hz/s，采样点数也相应地由58520变成21232（如下图所示）。为了得到重采样后的波形和频

5、谱，可运行如下语句：yy=fft(y);subplot(2,1,1);plot(y);title(重采样信号波形);subplot(2,1,2);plot(abs(yy); title(重采样信号频谱);运行后得到的波形、频谱图如下：（2）语音信号的频谱分析等间隔取8000个数据，首先画出语音信号的时域波形，然后对语音信号进行频谱分析。在MATLAB中，可以利用函数fft对信号进行快速傅里叶变换，得到信号的频谱特性(幅度和相位谱)。等间隔取16000个数据，画出语音信号的时域波形，然后对语音信号进行频谱分析。运用“四舍五入”法对y信号进行等间隔抽取，编写程序如下：m,n=size(y);z1=

6、zeros(1,8000);k1=1;for i=1:m/8000:m i=round(i); z1(k1)=y(i); k1=k1+1;endz2=zeros(1,16000);k2=1;for i=1:m/16000:m i=round(i); z2(k2)=y(i); k2=k2+1;end分别运行如下程序，输出各自对应波形、频谱图y1=fft(z1);subplot(3,2,3);plot(z1);title(等间隔采样8000点的波形);subplot(3,2,4);plot(abs(y1);title(等间隔采样8000点的频谱);y2=fft(z2);subplot(3,2,5)

7、;plot(z2);title(等间隔采样16000点的波形);subplot(3,2,6);plot(abs(y2);title(等间隔采样16000点的频谱);（3）设计FIR和IIR数字滤波器并画出频率响应根据有关的频谱特征，设计FIR和IIR数字滤波器。在MATLAB中，可以利用函数fir1设计FIR滤波器，利用函数butte、cheby1和ellip设计IIR滤波器；最后，利用MATLAB中的函数freqz画出各滤波器的频率响应。依题目要求，根据等间隔采样后频谱特征，分别设计出FIR低通滤波器、FIR高通滤波器、FIR带通滤波器、双线性滤波器等。 FIR低通滤波器用窗函数法设计FIR

8、低通滤波器，编写程序如下：fp=1000;fc=1200;as=100;ap=1;fs=22000;wp=2*fp/fs;wc=2*fc/fs;N=ceil(as-7.95)/(14.36*(wc-wp)/2)+1;beta=0.1102*(as-8.7);window=Kaiser(N+1,beta);b=fir1(N,wc,window);freqz(b,1,512,fs);用窗函数法设计的低通滤波器的幅度谱、相位谱 FIR高通滤波器用窗函数法设计FIR高通滤波器，编写程序如下：fs=22000;Wp=2*5000/fs;Ws=2*4800/fs;Ap=1;As=100;N=ceil(8*

9、pi/(Wp-Ws)+1;N=N+mod(N+1,2)+1;Wc=(Wp+Ws)/2/pi;b=fir1(N,Wc,high);omega=linspace(0,pi,512);freqz(b,1,omega);用窗函数法设计的高通滤波器的幅度谱、相位谱 FIR带通滤波器用窗函数法设计FIR带通滤波器，编写程序如下：fs=22000;Wp1=2*1200/fs;Wp2=2*3000/fs;Wc1=2*1000/fs;Wc2=2*3200/fs;Ap=1;As=100;W1=(Wp1+Wc1)/2;W2=(Wp2+Wc2)/2;wdth=min(Wp1-Wc1),(Wc2-Wp2);N=ceil

10、(11*pi/wdth)+1;b = fir1(N,W1 W2);freqz(b,1,512,fs);用窗函数法设计的带通滤波器的幅度谱、相位谱 双线性滤波器双线性滤波器的程序编写如下：fp=1000; fc=1200;As=100; Ap=1; fs=22050;wc=2*fc/fs;wp=2*fp/fs;n,wn=ellipord(wp,wc,Ap,As);b,a=ellip(n,Ap,As,wn);freqz(b,a,512,fs);用窗函数法设计的带通滤波器的幅度谱、相位谱（4）对信号进行滤波针对电话信道（最高3500Hz），对所有数据进行插值和抽取处理，把抽样率转变为7000Hz/s

11、，并进行频谱分析，得到幅度和相位谱。(此处必须把所有的数据参加FFT后进行低通滤波。可用（3）中的设计FIR或IIR滤波器进行时域滤波。在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。) FIR低通滤波器通过之前设计的FIR低通滤波器对z1、z2进行滤波并显示结果i1=fftfilt(b,z1);i2=fft(i1);i3=fftfilt(b,z2);i4=fft(i3);subplot(221); plot(i1); title(等间隔采用8000点波形经过滤波后信号的波形);subplot(222); plot(abs(i

12、2); title(等间隔采用8000点波形经过滤波后信号的频谱);subplot(223); plot(i3); title(等间隔采用16000点波形经过滤波后信号的波形);subplot(224); plot(abs(i4);title(等间隔采用16000点波形经过滤波后信号的频谱); FIR高通滤波器通过之前设计的FIR高通滤波器对z1、z2进行滤波并显示结果i1=fftfilt(b,z1);i2=fft(i1);i3=fftfilt(b,z2);i4=fft(i3);subplot(221); plot(i1);title(等间隔采用8000点波形经过滤波后信号的波形);subpl

13、ot(222); plot(abs(i2); title(等间隔采用8000点波形经过滤波后信号的频谱);subplot(223); plot(i3); title(等间隔采用16000点波形经过滤波后信号的波形);subplot(224); plot(abs(i4);title(等间隔采用16000点波形经过滤波后信号的频谱); FIR带通滤波器通过之前设计的FIR高通滤波器对z1、z2进行滤波并显示结果i1=fftfilt(b,z1);i2=fft(i1);i3=fftfilt(b,z2);i4=fft(i3);subplot(221); plot(i1); title(等间隔采用8000

14、点波形经过滤波后信号的波形);subplot(222); plot(abs(i2);title(等间隔采用8000点波形经过滤波后信号的频谱);subplot(223); plot(i3); title(等间隔采用16000点波形经过滤波后信号的波形);subplot(224); plot(abs(i4); title(等间隔采用16000点波形经过滤波后信号的频谱); 双线性滤波器通过之前设计的双线性滤波器对z1、z2进行滤波并显示结果i1=fftfilt(b,z1);i2=fft(i1);i3=fftfilt(b,z2);i4=fft(i3);subplot(221); plot(i1);

15、 title(等间隔采用8000点波形经过滤波后信号的波形);subplot(222); plot(abs(i2);title(等间隔采用8000点波形经过滤波后信号的频谱);subplot(223); plot(i3); title(等间隔采用16000点波形经过滤波后信号的波形);subplot(224); plot(abs(i4); title(等间隔采用16000点波形经过滤波后信号的频谱);（5）被污染语音信号的频谱分析把（4）处理后的所有数据储存为声音文件，与原始声音进行比较。（6）回放语音信号在MATLAB中，函数sound可以对声音进行回放。其调用格式为：sound(x,fs,

16、bits)可以感觉滤波前后的声音有变化。 根据以上各个滤波器的滤波结果，用sound函数可以听到声音，对比发现低通滤波器和带通滤波器的效果比较好。5.调试分析5.1等间隔采样程序的编写为了等间隔采样8000个数据，运行如下程序， x,fs,bits=wavread(d:5.wav);y=resample(x,8000,22050);m,n=size(y);a1=1:m/8000:m;b1=1:m/16000:m;a1=round(a1);b1=round(b1);a2=zeros(1,8000);b2=zeros(1,16000);for i=1:8000 a2(i)=y(i);endy1=f

17、ft(a2);subplot(2,2,1);plot(a2);subplot(2,2,2);plot(abs(y1);得到下图：与原信号对比发现等间隔采样8000点的波形全部集中在60008000之间，丢失了许多有用的信号，经过不断研究发现上述运行语句编写有误，通过不断地推敲修正，将程序改为：m,n=size(y);z1=zeros(1,8000);k1=1;for i=1:m/8000:m i=round(i); z1(k1)=y(i); k1=k1+1;end修改程序后即可得到正确的波形图：5.2 .m文件的命名把编写好的程序存为M文件（如下图所示），直接将该文件拉入到MATLAB软件运行

18、，运行结果出现错误（如下图所示）经过查找资料明白了matlab的.m文件保存的命名规则：1.文件名命名要用英文字符，第一个字符不能是数字2.文件名不要取为matlab的一个固有函数,m文件名的命名尽量不要是简单的英文单词，最好是由大小写英文/数字/下划线等组成。原因是简单的单词命名容易与matlab内部函数名同名，结果会出现一些莫名其妙的错误。3.文件存储路径一定为英文4.m文件起名不能为两个单词，如three phase，应该写成three_phase或者ThreePhase 对比自己M文件文件名可知是命名出现了错误，经过修改即可直接将M文件拉入MATLAB软件运行。5.3标点符号的使用在M

19、ATLAB软件中在中文输入法下如果输入标点符号就会显示如下错误：切换成英文输入法则可修正此错误。5.4滤波器参数的使用在使用IIR滤波器时，参数调整错误造成波形错误显示如下：6.结果分析与体会6.1总结 对采样频率、采样位数的认识经过不断地使用wavread、resample等函数，让我更深刻理解了采样位数与采样频率的区别： 采样位数（MATLAB中的bits）可以理解为声卡处理声音的解析度。这个数值越大，解析度就越高，录制和回放的声音就越真实。采样频率（MATLAB中的fs）是指录音设备在一秒钟内对声音信号的采样次数，采样频率越高声音的还原就越真实越自然。等间隔采样的方法等间隔取8000，1

20、6000个数据有许多方法，应该因地制宜，根据自己处理的语音信号的点数进行分析判断采用那种方法更适合。设待处理信号的点数为N， 若N比8000/16000的K（K为正整数）倍略多，则可从待处理信号的第一点开始采样，每隔K个点采样一次，共采样8000/16000个点，剩下的个别点数舍去；若N比8000/16000的K（K为正整数）倍略少，则可先在待处理信号的末尾补零，使信号的点数为8000/16000的K倍，然后从第一点开始采样，同样每隔K个点采样一次共采样8000/16000个点；若N比8000/16000的K（K为正整数）倍多很多，比8000/16000的K+1（K为正整数）倍少很多，如果用两

21、种方法势必丢失很多有用信号，这时可以用“四舍五入”法实现，从第一点开始采样，以（N/8000）的绝对值作为采样间隔，到第N点刚好是第8000个点，等间隔采样16000个点同理。因为离散信号不存在小数的n，所以当N/8000为小数时必须取绝对值，才能找到相对应的值。等间隔取样8000点和等间隔取样16000点之间的不同：对于y信号等间隔采样8000点和等间隔取样16000点直接造成了采样频率的不同，采样点数的多少也从一定程度上影响了音频信号的质量，用sound函数试听两种不同采样点数所形成的声音信号，对比可知：等间隔采样8000点信号的声音质量比等间隔采样16000点信号差，因为采样点数少使得原

22、信号的许多有用信号丢失，造成了语音信号的不完整。采样频率的选择采样频率是本次课程设计十分关键的一个量，贯穿始终，从第一步到最后一步都十分依赖于采样频率的取值。对于语音信号来说采样过程肯定会丢失一些数据,采样频率不够高还原后的信号就会失真，采样的点少了,如果没有达到信号的最高频率的两倍,那就有频谱混叠,恢复信号的时候和原信号就不一样了，因此在处理过程中要选择好采样频率.对滤波器的认识MATLAB提供了许多常见的滤波器，其中有很多常用的窗函数，包括hanning、hamming、blackman、Kaiser等几种，本次课程设计也是用窗函数法设计FIR滤波器，窗函数法设计FIR滤波器一般分为3个步

23、骤：第1步估计FIR滤波器阶数M（或长度N）。滤波器的阶数是根据技术指标确定的，可先设计滤波器,然后根据信号的频谱调整程序中的技术指标等参数。第2步确定所用的窗函数并计算出窗函数的值；第3步计算理想滤波器的单位脉冲响应并用窗函数将其截断即得所设计的FIR滤波器的hk。滤波器的设计需要不断地调试，通过对比不同技术指标滤波器处理下的声音不断地完善自己设计的滤波器。6.2不足本次课程设计的采样频率过小,直接导致大量的语音信号有效数据丢失（直接用wavread函数将连续信号采样后,采样频率足够大,不过经重采样后采样频率就小了），对比原始声音与经过重采样后的声音即可发现二者已有了很大的差别，这对之后的步

24、骤产生了很大的影响。6.3思考1从那几种数字滤波器的幅频特性曲线中可以观察到双线性变换法中和之间的非线性关系？ 2能否利用公式完成脉冲响应不变法的数字滤波器设计？6.4结束语本次课程设计培养了我们的学习兴趣，通过不断地学习、查资料懂得了不少课本上学不到的知识，其中少不了同学之间的互相帮助，老师的细心点评，集中体现了团结协作精神，我们投入热情、收获经验。通过这次的课程设计，巩固了以前C+的编程知识，强化了对数字信号处理、信号与系统、数据通信原理等课程原理的认识与把握，课本知识得到了升华。另外，通过本次课程设计，我学习了MATLAB软件信号处理方便的部分功能，深刻体会到该软件的强大，里面齐全的函数

25、以及配套完整的help语句帮助我们实现各种各样的要求。利用MATLAB的强大运算功能，基于MATLAB信号处理工具箱（Signal Processing Toolbox）的数字滤波器设计法可以快速有效的设计由软件组成的常规数字滤波器，设计方便、快捷，极大的减轻了工作量。在设计过程中可以对比滤波器特性，随时更改参数，以达到滤波器设计的最优化，以此来实现本次课程设计最重要的滤波环节。本次课程设计给了我诸多感想，脑子里折射出许多问题。像我们这种专业的学生确实应该多结合所学的理论知识，加强动手实践能力，多学习相关软件，多上网查询资料，巩固所学知识，课程设计中往往能引发我们的深思，发挥我们的想象力，促进

26、我们不断进步，给我们学习的动力。参考文献1高西全，丁玉美 编著 数字信号处理.西安电子科技大学出版社，20082陈后金，胡健，薛健 编著 信号与系统，高等教育出版社，20073刘波，文忠，曾涯 编著 MATLAB信号处理，电子工业出版社，20064万永革 编著，数字信号处理的MATLAB实现，科学出版社，20075李勇，徐震 编著 MATLAB辅助现代工程数字信号处理，西安电子科技大学出版社，20026张葛祥，李娜. 编著MATLAB仿真技术与应用 清华大学出版社，20037李海淘，邓樱 编著MATLAB程序设计教程 高等教育出版社 ,20028董长虹主编. 编著 MATLAB信号处理与应用

27、国防工业出版社，2005附录：x,fs,bits=wavread(d:5.wav);y=resample(x,8000,22050);m,n=size(y);z1=zeros(1,8000);k1=1;for i=1:m/8000:m i=round(i); z1(k1)=y(i); k1=k1+1;endz2=zeros(1,16000);k2=1;for i=1:m/16000:m i=round(i); z2(k2)=y(i); k2=k2+1;endsubplot(3,2,1); plot(y); title(原始信号波形);subplot(3,2,2); plot(abs(y); t

28、itle(原始信号频谱);y1=fft(z1);subplot(3,2,3);plot(z1);title(等间隔采样8000点的波形);subplot(3,2,4);plot(abs(y1);title(等间隔采样8000点的频谱);y2=fft(z2);subplot(3,2,5);plot(z2);title(等间隔采样16000点的波形);subplot(3,2,6);plot(abs(y2);title(等间隔采样16000点的频谱);fs=22000;Wp=2*5000/fs;Ws=2*4800/fs;Ap=1;As=100;N=ceil(8*pi/(Wp-Ws)+1;N=N+mo

29、d(N+1,2)+1;Wc=(Wp+Ws)/2/pi;b=fir1(N,Wc,high);omega=linspace(0,pi,512);freqz(b,1,omega);i1=fftfilt(b,z1);i2=fft(i1);i3=fftfilt(b,z2);i4=fft(i3);subplot(221); plot(i1); title(等间隔采用8000点波形经过滤波后信号的波形);subplot(222); plot(abs(i2); title(等间隔采用8000点波形经过滤波后信号的频谱);subplot(223); plot(i3); title(等间隔采用16000点波形经过滤波后信号的波形);subplot(224); plot(abs(i4); title(等间隔采用16000点波形经过滤波后信号的频谱);sound(i1,fs,bits);sound(i3,fs,bits);