数字信号课程设计-语音信号的处理与滤波.doc

广西科技大学（筹）数字信号处理课程设计设计题目：语音信号的处理与滤波设计过程：1.语音信号的采集利用电脑上的声卡和windows操作系统可以进行数字信号的采集。可以直接用电脑自带的录音工具进行录音。也可以运用matlab中的wavrecord()、wavwrite()、wavplay()、wavread()等函数进行语音的录制、播放、存储、读取。具体采集步骤如下： 用matlab工具编程录制一段时间为5秒的内容（例如“老师，新年快乐” ）的语音，存放在c:userspkldesktoptft.wav,并把存放录音文件的文件夹导入matlab的work目录中。录制程序如下： 程序： fs=16000;%设计抽样频率 channel=1;%设计录音通道 t=5; fprintf(按任意键开始%d秒录音:,t); pause;%暂停命令 fprintf(正在录音);y=wavrecord(t*fs,fs,channel,double);%录制语音 fprintf(录音结束n); wavwrite(y,fs,c:userspkldesktoptft.wav); fprintf(按任意键回放语音); pause;wavplay(y,fs);2. 语音信号的时频域分析 利用matlab中的“wavread”命令来读入已经录制好的语音信号，将它赋值给任意向量，再对其进行采样。然后可用matlab的绘图函数对时域和频域进行绘图分析。 运用下面程序进行采样和绘制时频域图，时域图为附件中的图1，频域图为附件中的图2。 程序： y,fs,nbits=wavread(tft.wav);%加载语音信号 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 原始信号频谱图);3. 采用窗函数法和双线性变换法设计滤波器所需设计的滤波器的性能指标为：（1） 低通滤波器性能指标：fb=1000hz,fc=1200hz,as=100db,ap=1db。（2） 高通滤波器性能指标：fb=5000hz,fc=4800hz,as=100db,ap=1db。（3） 带通滤波器性能指标：fp1=1200hz,fp2=3000hz,fsc1=1000hz,fsc2=3200hz,as=100db,ap=1db。(1） 用双线性变换法设计滤波器用双线性变换法设计滤波器就是先按低通滤波器的性能指标设计模拟低通滤波器，再转换成相应的过渡模拟滤波器。然后采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。这样可以避免产生频谱混叠现象。在双线性变换法中可用butter,cheby1和ellip设计相应iir滤波器。设计程序如下：低通滤波器设计程序： fp=1000; ft=5500;%取抽样频率 fs=1200; wp=2*pi*fp/ft;%通带频率 ws=2*pi*fs/ft;%截止频率 omgp=2*ft*tan(wp/2); omgs=2*ft*tan(ws/2); n,omgc=buttord(omgp,omgs,1,100,s); z,p,k=buttap(n); b=k*real(poly(z); a=real(poly(p); b,a=lp2lp(b,a,omgc); ba,aa=bilinear(b,a,ft); ha,w=freqz(ba,aa); plot(w/pi,abs(ha); title(图3 iir低通滤波器);附件图3为iir低通滤波器频谱图。 高通滤波器设计程序：fp=5000;ft=25000;%取抽样频率fs=4800;rp=1;as=100;wp=2*pi*fp/ft;%通带频率ws=2*pi*fs/ft;%截止频率omgp=2*ft*tan(wp/2);omgs=2*ft*tan(ws/2);n,omgc=ellipord(omgp,omgs,rp,as,s); %计算阶数n和截止频率z,p,k=ellipap(n,rp,as);b=k*real(poly(z);a=real(poly(p);b,a=lp2hp(b,a,omgc);ba,aa=bilinear(b,a,ft);ha,w=freqz(ba,aa);plot(w/pi,abs(ha);title(图4 iir高通滤波器);附件图4为iir高通滤波器频谱图。 带通滤波器设计程序：fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;as=100;rp=1;wp1=2*pi*fp1/ft;wp2=2*pi*fp2/ft;ws1=2*pi*fs1/ft;ws2=2*pi*fs2/ft;omgp1=2*ft*tan(wp1/2);omgp2=2*ft*tan(wp2/2);omgp=omgp1,omgp2;omgs1=2*ft*tan(ws1/2);omgs2=2*ft*tan(ws2/2);omgs=omgs1,omgs2;bw=omgp2-omgp1;w0=sqrt(omgp1*omgp2);%通带宽和中心频率n,omgn=cheb1ord(omgp,omgs,rp,as,s);z,p,k=cheb1ap(n,rp);b=k*real(poly(z);a=real(poly(p);b,a=lp2bp(b,a,w0,bw);ba,aa=bilinear(b,a,ft);ha,w=freqz(ba,aa);plot(w/pi,abs(ha);title(图5iir带通滤波器);附件图5为iir高通滤波器频谱图。 (2） 用窗函数法设计滤波器 窗函数设计方法主要有矩形窗、巴特列特窗、汉宁窗、海明窗、布拉克曼窗凯泽窗等。选择合适窗设计fir滤波器，这里按照设置指标选择凯泽窗和布莱克曼窗进行设计。fir低通滤波器设计程序：ft=5500;fpts=1000 1200;dev=0.01 0.05;%通带与阻带波动mag=1 0;n21,wn21,beta,ftype=kaiserord(fpts,mag,dev,ft);b21=fir1(n21,wn21,kaiser(n21+1,beta);h,w=freqz(b21,1);plot(w/pi,abs(h);title(图6 fir低通滤波器);频率响应图为：图6fir高通滤波器设计程序：ft=25000;fpts=4800 5000;dev=0.01 0.05;mag=1 0;n22,wn22,beta,ftype=kaiserord(fpts,mag,dev,ft);b22=fir1(n22,wn22,high,kaiser(n22+1,beta);h,w=freqz(b22,1);plot(w/pi,abs(h);title(图7 fi高通滤波器);频率响应图为：图7fir带通滤波器的设计程序：fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;ws1=fs1/(ft/2)*pi;ws2=fs2/(ft/2)*pi;wp1=fp1/(ft/2)*pi;wp2=fp2/(ft/2)*pi;wp=wp1,wp2;ws=ws1,ws2;deltaw=wp1-ws1;n0=ceil(11*pi/deltaw);n=n0+mod(n0+1,2);windows=(blackman(n);wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=ideallp(wc2,n)-ideallp(wc1,n);b=hd.*windows;h,w=freqz(b,1);plot(w/pi,abs(h); title(图8 fir带通滤波器); 频率响应图为：图84. 对“tft.wav”语音信号进行滤波 分别采用设计好的低通，高通，带通滤波器对“tft.wav”语音信号进行滤波，用matlab画出滤波后的时频域图，并回放录播后的语音信号。（1） fir低通滤波程序：ft=5500;fpts=1000 1200;dev=0.01 0.05;mag=1 0;n21,wn21,beta,ftype=kaiserord(fpts,mag,dev,ft);b21=fir1(n21,wn21,kaiser(n21+1,beta);h,w=freqz(b21,1);y,fs,nbits=wavread(tft.wav);%加载语音信号f=filter(b21,1,y);wavplay(f);f1=fft(f);subplot(2,1,1)plot(f);title(低通滤波后的波形);subplot(2,1,2)plot(abs(f1);title(低通滤波后的频谱图)运行结果图为附件：图9（2） fir高通滤波程序：ft=25000;fpts=4800 5000;dev=0.01 0.05;mag=1 0;n22,wn22,beta,ftype=kaiserord(fpts,mag,dev,ft);b22=fir1(n22,wn22,high,kaiser(n22+1,beta);h,w=freqz(b22,1);y,fs,nbits=wavread(tft.wav);%加载语音信号f=filter(b22,1,y);wavplay(f);f1=fft(f);subplot(2,1,1)plot(f);title(高通滤波后的波形);subplot(2,1,2)plot(abs(f1);title(高通滤波后的频谱图)运行结果图为附件：图10（3） fir带通滤波程序：fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;ws1=fs1/(ft/2)*pi;ws2=fs2/(ft/2)*pi;wp1=fp1/(ft/2)*pi;wp2=fp2/(ft/2)*pi;wp=wp1,wp2;ws=ws1,ws2;deltaw=wp1-ws1;n0=ceil(11*pi/deltaw);n=n0+mod(n0+1,2);windows=(blackman(n);wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=ideallp(wc2,n)-ideallp(wc1,n);b=hd.*windows;h,w=freqz(b,1); y,fs,nbits=wavread(tft.wav);%加载语音信号 f=filter(b,1,y); wavplay(f); f1=fft(f); subplot(2,1,1) plot(f); title(带通滤波后的波形); subplot(2,1,2) plot(abs(f1); title(带通滤波后的频谱图) 运行结果图为附件：图11（4） iir低通滤波程序：fp=1000; ft=5500;%取抽样频率 fs=1200; wp=2*pi*fp/ft;%通带频率 ws=2*pi*fs/ft;%截止频率 omgp=2*ft*tan(wp/2); omgs=2*ft*tan(ws/2); n,omgc=buttord(omgp,omgs,1,100,s); z,p,k=buttap(n); b=k*real(poly(z); a=real(poly(p); b,a=lp2lp(b,a,omgc); ba,aa=bilinear(b,a,ft);ha,w=freqz(ba,aa);y,fs,nbits=wavread(tft.wav);%加载语音信号f=filter(ba,aa,y);wavplay(f);f1=fft(f);subplot(2,1,1)plot(f);title(低通滤波后的波形);subplot(2,1,2)plot(abs(f1);title(低通滤波后的频谱图)运行结果图为附件：图91（5） iir高通滤波程序： fp=5000; ft=25000;%取抽样频率 fs=4800; rp=1;as=100; wp=2*pi*fp/ft;%通带频率 ws=2*pi*fs/ft;%截止频率 omgp=2*ft*tan(wp/2); omgs=2*ft*tan(ws/2); n,omgc=ellipord(omgp,omgs,rp,as,s); %计算阶数n和截止频率 z,p,k=ellipap(n,rp,as); b=k*real(poly(z); a=real(poly(p); b,a=lp2hp(b,a,omgc); ba,aa=bilinear(b,a,ft);ha,w=freqz(ba,aa);y,fs,nbits=wavread(tft.wav);%加载语音信号f=filter(ba,aa,y);wavplay(f);f1=fft(f);subplot(2,1,1)plot(f);title(iir高通滤波后的波形);subplot(2,1,2)plot(abs(f1);title(iir高通滤波后的频谱图)运行结果图为附件：图101（6） iir带通滤波程序：fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;as=100;rp=1;wp1=2*pi*fp1/ft;wp2=2*pi*fp2/ft;ws1=2*pi*fs1/ft;ws2=2*pi*fs2/ft;omgp1=2*ft*tan(wp1/2);omgp2=2*ft*tan(wp2/2);omgp=omgp1,omgp2;omgs1=2*ft*tan(ws1/2);omgs2=2*ft*tan(ws2/2);omgs=omgs1,omgs2;bw=omgp2-omgp1;w0=sqrt(omgp1*omgp2);%通带宽和中心频率n,omgn=cheb1ord(omgp,omgs,rp,as,s);z,p,k=cheb1ap(n,rp);b=k*real(poly(z);a=real(poly(p);b,a=lp2bp(b,a,w0,bw);ba,aa=bilinear(b,a,ft);ha,w=freqz(ba,aa);y,fs,nbits=wavread(tft.wav);%加载语音信号f=filter(ba,aa,y);wavplay(f);f1=fft(f);subplot(2,1,1)plot(f);title(iir带通滤波后的波形);subplot(2,1,2)plot(abs(f1);title(iir带通滤波后的频谱图)运行结果图为附件：图1115. 异性语音采集分析 运用相关matlab函数和所学的相关知识，采集同样内容的异性的语音信号，对它的时频域进行分析处理，比较两种相同内容而不同性别语音信号的特点。步骤如下：（1） 异性语音录制运用matlab中的wavrecord()、wavwrite()、wavplay()等函数进行语音的录制、播放、存储。把录制好的语音文件存为“c:userspkldesktoptft2.wav”，具体程序如下： fs=16000;%设计抽样频率 channel=1;%设计录音通道 t=5; fprintf(按任意键开始%d秒录音:,t); pause;%暂停命令 fprintf(正在录音); y=wavrecord(t*fs,fs,channel,double);%录制语音 fprintf(录音结束n); wavwrite(y,fs,c:userspkldesktoptft2.wav); fprintf(按任意键回放语音); pause; wavplay(y,fs);（2） 异性语音的采样分析 运用matlab中的wavread()对录制好的异性语音信号进行加载采样，然后用matlab的相关函数对采集的语音信号进行时频域画图分析，作好的图为图12。具体程序如下： y,fs,nbits=wavread(tft2.wav);%加载语音信号 wavplay(y,fs);%回放语音信号 n=length(y);%语音信号长度 y=fft(y,n);%用傅里叶变换求频域 subplot(2,1,1); plot(y); title(异性语音原始信号时域波形图); subplot(2,1,2); plot(abs(y); title(异性语音原始信号频谱图);6. 语音信号加噪、滤噪处理分析运用matlab相关函数录制一段相同长度的背景噪声，再把此噪声叠加到原始语音信号中，然后再选择合适的滤波器对加噪语音信号进行滤波。步骤如下：（1） 噪声录制与叠加运用matlab随即录制一段与语音信号相同长度的噪声，再叠加到原始语音信号上，把加噪后的语音信号存为“c:userspkldesktoptft3.wav”。实现程序如下：y1,fs,nbits=wavread(tft.wav);y1=fft(y1);n=length(y1);noise=0.01*randn(n,1);%随机加同样长度的噪声y2=y1+noise;y2=fft(y2);sound(y2)wavwrite(y2,fs,c:userspkldesktoptft3.wav);subplot(2,2,1);plot(y1);title(叠加前语音信号波形);subplot(2,2,3);plot(abs(y1);title(叠加前语音信号频谱图);subplot(2,2,2);plot(y2);title(叠加后语音信号波形);subplot(2,2,4);plot(abs(y2);title(叠加后语音信号频谱图);程序所得图形为附件：图13（2） 滤除噪声通过（1）中的加噪前后时频域的分析，滤除此噪声应用低通滤波器进行滤波，实现程序如下：ft=5500;fpts=1000 1200;dev=0.01 0.05;%通带与阻带波动mag=1 0;n21,wn21,beta,ftype=kaiserord(fpts,mag,dev,ft);b21=fir1(n21,wn21,kaiser(n21+1,beta);h,w=freqz(b21,1);y,fs,nbits=wavread(tft3.wav);%加载语音信号y=fft(y);f=filter(b21,1,y);f=fft(f);subplot(2,2,1)plot(y);title(低通噪声滤波前的波形);subplot(2,2,3)plot(abs(y);title(低通噪声滤波前的频谱图)subplot(2,2,2)plot(f);title(低通噪声滤波后的波形);subplot(2,2,4)plot(abs(f);title(低通噪声滤波后的频谱图)运行结果图为附件：图147. 信号处理系统界面设计建立gui文件：file-new-gui:如图15，后选择blan gui(default),得到如图16的gui开发环境。布置好各控件，为控件编写程序来实现两数相加的功能。控件设计如图17,各按钮设计程序如下：gui设计程序：function varargout = untitled1(varargin)gui_singleton = 1;gui_state = struct(gui_name, mfilename, . gui_singleton, gui_singleton, . gui_openingfcn, untitled1_openingfcn, . gui_outputfcn, untitled1_outputfcn, . gui_layoutfcn, , . gui_callback, );if nargin & ischar(varargin1) gui_state.gui_callback = str2func(varargin1);endif nargout varargout1:nargout = gui_mainfcn(gui_state, varargin:);else gui_mainfcn(gui_state, varargin:);endfunction untitled1_openingfcn(hobject, eventdata, handles, varargin)handles.output = hobject; guidata(hobject, handles);function varargout = untitled1_outputfcn(hobject, eventdata, handles) varargout1 = handles.output;function pushbutton1_callback(hobject, eventdata, handles)global signal;sound(signal)function pushbutton2_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global y;n=length(signal);y=fft(signal,n);plot(signal);title(原始信号波形);axes(handles.axes2)plot(abs(y);title(原始信号频谱);function pushbutton3_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global si;n=length(signal);noise=0.01*randn(n,1);si=signal+noise;sound(si)function pushbutton4_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global si;plot(si);title(加噪语音信号的时域波形);axes(handles.axes2)cla;s=fft(si);plot(abs(s);title(加噪语音信号的频域波形);function pushbutton5_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global b21;ft=5500;fpts=1000 1200;dev=0.01 0.05;mag=1 0;n21,wn21,beta,ftype=kaiserord(fpts,mag,dev,ft);b21=fir1(n21,wn21,kaiser(n21+1,beta);h,w=freqz(b21,1);plot(w/pi,abs(h);title(fir 低通滤波器);function pushbutton6_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global b;fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;ws1=fs1/(ft/2)*pi;ws2=fs2/(ft/2)*pi;wp1=fp1/(ft/2)*pi;wp2=fp2/(ft/2)*pi;wp=wp1,wp2;ws=ws1,ws2;deltaw=wp1-ws1;n0=ceil(11*pi/deltaw);n=n0+mod(n0+1,2);windows=(blackman(n);wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=ideallp(wc2,n)-ideallp(wc1,n);b=hd.*windows;h,w=freqz(b,1);plot(w/pi,abs(h);title(fir 带通滤波器);function pushbutton7_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global b22;ft=25000;fpts=4800 5000;dev=0.01 0.05;mag=1 0;n22,wn22,beta,ftype=kaiserord(fpts,mag,dev,ft);b22=fir1(n22,wn22,high,kaiser(n22+1,beta);h,w=freqz(b22,1);plot(w/pi,abs(h);title(fir 高通滤波器);function pushbutton8_callback(hobject, eventdata, handles)axes(handles.axes2)cla;global b1;global ba1;global aa1;global omgp;global omgs;fp=1000;ft=5500;fs=1200;wp=2*pi*fp/ft;ws=2*pi*fs/ft;omgp=2*ft*tan(wp/2);omgs=2*ft*tan(ws/2);n,omgc=buttord(omgp,omgs,1,100,s); z,p,k=buttap(n);b=k*real(poly(z);a=real(poly(p);b1,a1=lp2lp(b,a,omgc);ba1,aa1=bilinear(b1,a1,ft);ha,w=freqz(ba1,aa1);plot(w/pi,abs(ha);title(iir 低通滤波器);function pushbutton9_callback(hobject, eventdata, handles)axes(handles.axes2)cla;global b2;global ba2;global aa2;fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;as=100;rp=1;wp1=2*pi*fp1/ft;wp2=2*pi*fp2/ft;ws1=2*pi*fs1/ft;ws2=2*pi*fs2/ft;omgp1=2*ft*tan(wp1/2);omgp2=2*ft*tan(wp2/2);omgp=omgp1,omgp2;omgs1=2*ft*tan(ws1/2);omgs2=2*ft*tan(ws2/2);omgs=omgs1,omgs2;bw=omgp2-omgp1;w0=sqrt(omgp1*omgp2);%通带宽和中心频率n,omgn=cheb1ord(omgp,omgs,rp,as,s);z,p,k=cheb1ap(n,rp);b2=k*real(poly(z);a2=real(poly(p);b2,a2=lp2bp(b2,a2,w0,bw);ba2,aa2=bilinear(b2,a2,ft);ha,w=freqz(ba2,aa2);plot(w/pi,abs(ha);title(iir 带通滤波器);function pushbutton10_callback(hobject, eventdata, handles)axes(handles.axes2)cla;global b3;global ba3;global aa3;global omgp;global omgs;fp=5000;ft=25000;%取抽样频率fs=4800;rp=1;as=100;wp=2*pi*fp/ft;%通带频率ws=2*pi*fs/ft;%截止频率omgp=2*ft*tan(wp/2);omgs=2*ft*tan(ws/2);n,omgc=ellipord(omgp,omgs,rp,as,s); %计算阶数n和截止频率z,p,k=ellipap(n,rp,as);b3=k*real(poly(z);a3=real(poly(p);b3,a3=lp2hp(b3,a3,omgc);ba3,aa3=bilinear(b3,a3,ft);ha,w=freqz(ba3,aa3);plot(w/pi,abs(ha);title(iir高通滤波器);function pushbutton11_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global b21;global si;z21=fftfilt(b21,si);sound(z21)m21=fft(z21);plot(z21);title(滤波后的信号波形);axes(handles.axes2)cla;plot(abs(m21),r);title(滤波后信号频谱);function pushbutton12_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global b;global si;z22=fftfilt(b,si);sound(z22)m22=fft(z22);plot(z22);title(滤波后的信号波形);axes(handles.axes2)cla;plot(abs(m22),r);title(滤波后信号频谱);function pushbutton13_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global b22;global si;z23=fftfilt(b22,si);sound(z23)m23=fft(z23);plot(z23);title(滤波后的信号波形);axes(handles.axes2)cla;plot(abs(m23),r);title(滤波后信号频谱);function pushbutton14_callback(hobject, eventdata, handles)axes(handles.axes1)cla;global signal;global b1;global ba1;global aa1;global si;z11=filter(ba1,aa1,si);sound(z11);m11=fft(z11);plot(z11);title(滤波后的信号波形);axes(handles.axes2)cla;plot(abs(m11),r);title(滤波后信号频谱);function pushbutton15_callback(hobje