应用Matlab对语音信号进行频谱分析及滤波.doc

数字信号处理课程设计报告书应用Matlab对语音信号进行频谱分析及滤波题目：应用Matlab对语音信号进行频谱分析及滤波一、课程设计内容录制一段个人自己的语音信号，画出此语音信号的时域波形和频谱图；给定滤波器的性能指标，设计IIR和FIR滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，设计一个信号处理系统界面。二、课程设计目的1、进一步巩固数字信号处理的基本概念、理论、分析方法和实现方法；2、增强应用Matlab语言编写数字信号处理的应用程序及分析、解决实际问题的能力。三、设计内容及步骤1、语音信号的采集利用matlab的内置函数，设计一个录音器，录制一段自己的话音，时间在1s内。然后在Matlab软件平台下，利用函数wavread提取录制的语音信号，记住采样频率和采样点数。说明：由于录音程序录制后用wavwrite和wavrecord时一直出现错误，无法运行，在这里用filename,filepath=uigetfile(.wav,Open wav file);y,fs,nbits=wavread(filepath,filename);格式打开本地一个已录制好的wav文件进行分析。2 、语音信号的频谱分析首先画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性，从而加深对频谱特性的理解。3、设计数字滤波器和画出其频率响应用窗函数法和双线性变换法设计下面要求的3种滤波器。 (1)低通滤波器性能指标通带截止频率1000Hz，阻带截止频率1200 Hz，通带最大衰减1dB，阻带最小衰减=100dB。(2)高通滤波器性能指标通带截止频率5000Hz，阻带截止频率3000 通带衰减1dB，阻带衰减100dB。(3)带通滤波器性能指标，通带范围是（1200Hz，3000Hz），在此两频率处衰减不大于1dB，在1000Hz和3200Hz频率处衰减不小于100dB。在Matlab中，可以利用函数fir1设计FIR滤波器，可以利用函数butte，cheby1和ellip设计IIR滤波器；利用Matlab中的函数freqz画出各滤波器的频率响应。4 、用滤波器对信号进行滤波用自己设计的各滤波器分别对采集的信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。5、 比较滤波前后语音信号的波形及频谱在一个窗口同时画出滤波前后的波形及频谱。6 、回放语音信号在Matlab中，函数sound可以对声音进行回放。其调用格式：sound(x，fs，bits)；可以感觉滤波前后的声音有变化。7 、设计系统界面本系统界面主要采用ifelse语句和menu菜单函数进行设计，如下图： 四、结果分析1、原始语音信号采样后的时域波形及FFT变换后频谱2、FIR滤波器及原始信号通过FIR后的时域波形频谱变化(1)FIR低通滤波器FIR低通滤波器频率响应图：分析：由图可以看出，该滤波器过滤掉了高频部分。(2)FIR高通滤波器FIR高通滤波器频率响应图：分析：由图可以看出，该滤波器将低频成分滤掉。(3)FIR带通滤波器FIR带通滤波器频率响应图：分析：由上图看出该滤波器将低频和高频成分滤掉，剩下通带里的频率。3、IIR滤波器及原始信号通过FIR后的时域波形频谱变化(1)IIR低通滤波器分析：题目要求,最大衰减AP=1dB=0.89时fb=1000Hz，由图可看出，基本符合，而最小衰减为As=100dB= 0.00005时fc=1200Hz，基本符合要求。分析：由图可以看出，该滤波器过滤掉了高频部分。（2）IIR高通滤波器分析：题目要求,最大衰减AP=1dB=0.89时fb=5000Hz，由图可看出，基本符合，而最小衰减为As=100dB= 0.00005时fc=4800Hz，基本符合要求。分析：由图可以看出，该滤波器将低频成分滤掉。（2）IIR带通滤波器IIR带通滤波器频率响应图：分析：题目要求,最大衰减AP=1dB=0.89时fb1=1200Hz，fb2=3000Hz由图可看出，基本符合，而最小衰减为As=100dB= 0.00005时fc1=1000Hz，fc2=3200Hz基本符合。分析：上图很明显看出该滤波器将低频和高频成分滤掉了，只剩下通带里的频率。五、心得体会第一次做有关于Matlab的课程设计，虽然在以前也有类似的课题，但感觉上我对数字信号这门课程掌握的远远不够。首先是Matlab的知识不扎实，其次是滤波器的设计让我觉得有些困难。题目要求是要学生自己写录音程序，但在设计程序过程中，说明：由于录音程序录制后调用wavwrite和wavrecord时一直出现错误，无法运行，花费了不少时间，无奈之下只好用一下程序代替：filename,filepath=uigetfile(.wav,Open wav file);y,fs,nbits=wavread(filepath,filename);附上一段写过的录音程序：fs=16000; %取样频率t=1; %录音时间fprintf(Press any key to start %g seconds of recording.n,duration);pause;fprintf(Recording.n);y=wavrecord(t*fs,fs); %duration*fs 是总的采样点数fprintf(Finished recording.n);fprintf(Press any key to play the recording.n);pause;wavplay(y,fs);由于时间仓促，这次课程设计在系统界面设计上采用Menu以及两种滤波器的六种基本形式，比较遗憾。尽管如此，通过做应用Matlab对语音信号进行谱分析及滤波这个课程设计，加深了我对数字信号处理理论知识，尤其是采样频率，频谱特性和数字滤波器的设计；另一方面，使我在原来熟悉控制系统在Matlab运用的基础上进一步学习Matlab软件在信号分析和处理的运用。分析课题题目：应用Matlab对语音信号进行谱分析及滤波，考虑到几个方面的内容：（1）加强对Matlab基础知识的学习，尤其是Matlab在信号处理方面的学习，这就需要多方面查找资料，如上互联网，图书馆等；（2）加强学习信号处理知识理论，尤其是FFT频谱变换理论和各种数字滤波器设计理论；（3）要突出重点，其重点是谱分析和数字滤波器的设计和分析结果。 六、设计程序clear ;clc;filename,filepath=uigetfile(.wav,Open wav file);y,fs,nbits=wavread(filepath,filename); %从电脑文件夹选择wav音频文件b=menu(请选择选项,原始信号采样后时域图和频谱图,FIR滤波器,IIR滤波器,退出);while(b=4)if b=1 temp=menu(请选择选项,播放原始语音,原始语音时域图,原始语音做FFT变换后频谱图,返回); if temp=1 %播放语音 sound(y); elseif temp=2 %画出原始语音时域图 figure(11) plot(y); xlabel(时间（ms); ylabel(幅值); title(原始语音时域图); elseif temp=3 %画出原始语音做FFT变换后频谱图 Y1=fft(y); Y=abs(Y1); figure(12) plot(Y); xlabel(频率); ylabel(幅值); title(原始语音做FFT变换后频谱图); else endend if b=2 temp=menu(请选择选项,FIR低通滤波器,FIR高通滤波器,FIR带通滤波器,返回); if temp=1 %FIR低通滤波器 fs=10000; wp=2*pi*1000/fs; wst=2*pi*1200/fs; Rp=1; Rs=100; wdelta=wst-wp; N=ceil(8*pi/wdelta); %取整 wn=(wp+wst)/2; b,a=fir1(N,wn/pi,hamming(N+1); %选择窗函数，并归一化截止频率 figure(21) freqz(b,a,512); title(FIR低通滤波器); y1=filter(b,a,y); figure(22) subplot(2,1,1) plot(y) title(FIR低通滤波器滤波前的时域波形); xlabel(时间（ms); ylabel(幅值); subplot(2,1,2) plot(y1); title(FIR低通滤波器滤波后的时域波形); xlabel(时间（ms); ylabel(幅值); sound(y1,8000); %播放滤波后的语音信号 F0=fft(y1,1024); f=fs*(0:511)/1024; figure(23) y2=fft(y,1024); subplot(2,1,1); plot(f,abs(y2(1:512); title(FIR带通滤波器滤波前的频谱) xlabel(频率/Hz); ylabel(幅值); subplot(2,1,2) F2=plot(f,abs(F0(1:512); title(FIR低通滤波器滤波后的频谱); xlabel(频率/Hz); ylabel(幅值); elseif temp=2 %FIR高通滤波器 fs=22050; wp=2*pi*5000/fs; wst=2*pi*3000/fs; Rp=1; Rs=100; wdelta=wp-wst; N=ceil(8*pi/wdelta); %取整 wn=(wp+wst)/2; b,a=fir1(N,wn/pi,high); figure(24) freqz(b,a,512); title(FIR高通滤波器); y1=filter(b,a,y); figure(25) subplot(2,1,1) plot(y) title(FIR高通滤波器滤波前的时域波形); xlabel(时间（ms); ylabel(幅值); subplot(2,1,2) plot(y1); title(FIR高通滤波器滤波后的时域波形); xlabel(时间（ms); ylabel(幅值); sound(y1,8000); %播放滤波后的语音信号 F0=fft(y1,3000); f=fs*(0:511)/1024; figure(26) y2=fft(y,1024); subplot(2,1,1); plot(f,abs(y2(1:512); title(FIR高通滤波器滤波前的频谱) xlabel(频率/Hz); ylabel(幅值); subplot(2,1,2) plot(f,abs(F0(1:512); title(FIR高通滤波器滤波后的频谱) xlabel(频率/Hz); ylabel(幅值); elseif temp=3 %FIR带通滤波器 fs=22050; wp1=2*pi*1200/fs;wp2=2*pi*3000/fs; ws1=2*pi*1000/fs;ws2=2*pi*3200/fs; Rp=1; Rs=100; wp=(wp1+ws1)/2;ws=(wp2+ws2)/2; wdelta=wp1-ws1; N=ceil(8*pi/wdelta); %取整 wn=wp ws; b,a=fir1(N,wn/pi,bandpass); figure(27) freqz(b,a,512); title(FIR带通滤波器); y1=filter(b,a,y); figure(28) subplot(2,1,1) plot(y); xlabel(时间（ms); ylabel(幅值); title(FIR带通滤波器滤波前的时域波形); subplot(2,1,2) plot(y1); title(FIR带通滤波器滤波后的时域波形); xlabel(时间（ms); ylabel(幅值); sound(y1,8000); %播放滤波后的语音信号 F0=fft(y1,3000); f=fs*(0:511)/1024; figure(29) y2=fft(y,1024); subplot(2,1,1); plot(f,abs(y2(1:512); title(FIR带通滤波器滤波前的频谱) xlabel(频率/Hz); ylabel(幅值); subplot(2,1,2) plot(f,abs(F0(1:512); title(FIR带通滤波器滤波后的频谱) xlabel(频率/Hz); ylabel(幅值); else endend if b=3 temp=menu(请选择选项,IIR低通滤波器,IIR高通滤波器,IIR带通滤波器,返回); if temp=1 %IIR低通滤波器 fs=22050; Ts=1/fs; R1=10; wp=2*pi*1000/fs; ws=2*pi*1200/fs; Rp=1; Rl=100; wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标 ws1=2/Ts*tan(ws/2); N,Wn=buttord(wp1,ws1,Rp,R1,s); %选择滤波器的最小阶数 Z,P,K=buttap(N); %创建butterworth模拟滤波器 Bap,Aap=zp2tf(Z,P,K); b,a=lp2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 H,W=freqz(bz,az); figure(31) plot(W*fs/(2*pi),abs(H) grid xlabel(频率Hz); ylabel(频率响应幅度); title(IIR低通滤波器); f1=filter(bz,az,y); figure(32) subplot(2,1,1) %画出滤波前的时域图 plot(y); xlabel(时间（ms); ylabel(幅值); title(IIR低通滤波器滤波前的时域波形); subplot(2,1,2) plot(f1); %画出滤波后的时域图 title(IIR低通滤波器滤波后的时域波形); xlabel(时间（ms); ylabel(幅值); sound(f1,9000); %播放滤波后的信号 F0=fft(f1,3000); f=fs*(0:511)/1024; figure(33) y2=fft(y,1024); subplot(2,1,1); plot(f,abs(y2(1:512); %画出滤波前的频谱图 title(IIR低通滤波器滤波前的频谱) xlabel(频率/Hz); ylabel(幅值); subplot(2,1,2) f=fs*(0:511)/1024; F1=plot(f,abs(F0(1:512); %画出滤波后的频谱图 title(IIR低通滤波器滤波后的频谱) xlabel(频率/Hz); ylabel(幅值); elseif temp=2 %IIR高通滤波器 fs=22050; Fs=22050; Ts=1/Fs; R1=50; Wp=2*pi*5000/fs; Ws=2*pi*3000/fs; Rp=1; Rl=100; Wp1=2/Ts*tan(Wp/2); %将模拟指标转换成数字指标 Ws1=2/Ts*tan(Ws/2); N,Wn=cheb2ord(Wp1,Ws1,Rp,Rl,s); %选择滤波器的最小阶数 Z,P,K=cheb2ap(N,Rl); %创建切比雪夫模拟滤波器 Bap,Aap=zp2tf(Z,P,K); b,a=lp2hp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 H,W=freqz(bz,az); %绘制频率响应曲线 figure(34) plot(W*fs/(2*pi),abs(H); grid xlabel(频率Hz); ylabel(频率响应幅度); title(IIR高通滤波器); f1=filter(bz,az,y); figure(35) subplot(2,1,1) plot(y); %画出滤波前的时域图 xlabel(时间（ms); ylabel(幅值); title(IIR高通滤波器滤波前的时域波形); subplot(2,1,2) xlabel(时间（ms); ylabel(幅值); plot(f1); xlabel(时间（ms); ylabel(幅值); title(IIR高通滤波器滤波后的时域波形); sound(f1,8000); %播放滤波后的信号 F0=fft(f1,1024); figure(36) y2=fft(y,1024); subplot(2,1,1); plot(f,abs(y2(1:512); %画出滤波前的频谱图 title(IIR高通滤波器滤波前的频谱) xlabel(频率/Hz); ylabel(幅值); subplot(2,1,2) f=fs*(0:511)/1024; plot(f,abs(F0(1:512); %画出滤波后的频谱图 title(IIR高通滤波器滤波后的频谱) xlabel(频率/Hz); ylabel(幅值); elseif temp=3 %IIR带通滤波器 Fs=22050; Ts=1/Fs; R1=30; fb1=1200; fb2=3000; fc1=1000; fc2=3200; fs=22050; W1=2*fb1*pi/fs; W2=2*fc1*pi/fs; W3=2*fb2*pi/fs; W4=2*fc2*pi/fs; Wp=W1,W3; Ws=W2,W4; Rp=1; Rl=100; Wp1=2/Ts*tan(Wp/2); %将模拟指标转换成数字指标 Ws1=2/Ts*tan(Ws/2); N,Wn=cheb2ord(Wp1,Ws1,Rp,R1,s); %选择滤波器的最小阶数 Z,P