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

课程设计报告 课程名称： 数字信号处理 课题名称： 语音信号的处理与滤波 姓 名： 学 号： 院 系： 专业班级： 指导教师： 完成日期： 2013年7月2日 目录第1部分 课程设计报告3一设计目的3二设计内容3三设计原理3四具体实现51. 录制一段声音52. 巴特沃斯滤波器的设计83. 将声音信号送入滤波器滤波134. 语音信号的回放195. 男女语音信号的频谱分析196. 噪声的叠加和滤除22五 结果分析27第2部分 课程设计总结28一 参考文献28第1部分 课程设计报告一设计目的综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。二设计内容录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除； 三设计原理1.在Matlab软件平台下，利用函数wavrecord()，wavwrite(),wavread(),wavplay()对语音信号进行录制，存储，读取，回放。2.用y=fft(x)对采集的信号做快速傅立叶变换，并用h1,w=freqz(h)进行DTFT变换。3.掌握FIR DF线性相位的概念，即线性相位对及零点的约束，了解四种FIR DF的频响特点。4.在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波。5.抽样定理 连续信号经理想抽样后时域、频域发生的变化（理想抽样信号与连续信号频谱之间的关系） 理想抽样信号能否代表原始信号、如何不失真地还原信号即由离散信号恢复连续信号的条件（抽样定理） 理想采样过程描述：时域描述：频域描述：利用傅氏变换的性质，时域相乘频域卷积，若 则有与的关系：理想抽样信号的频谱是连续信号频谱的周期延拓,重复周期为Ws(采样角频率)。如果：即连续信号是带限的，且信号最高频率不超过抽样频率的二分之一，则可不失真恢复。奈奎斯特采样定理：要使实信号采样后能够不失真还原，采样频率必须大于信号最高频率的两倍：四具体实现1.录制一段声音1.1录制并分析在MATLAB中用wavrecord、wavread、wavplay、wavwrite对声音进行录制、读取、回放、存储。程序如下：Fs=8000; %抽样频率time=3; %录音时间fprintf(按Enter键录音%ds,time); %文字提示pause; %暂停命令fprintf(录音中.);x=wavrecord(time*Fs,Fs,double); %录制语音信号fprintf(录音结束); %文字提示fprintf(按Enter键回放录音); pause; %暂停命令wavplay(x,Fs); %按任意键播放语音信号wavwrite(x,Fs,C:UsersacerDesktop数字信号sound.wav); %存储语音信号N=length(x);%返回采样点数df=fs/N;%采样间隔n1=1:N/2;f=(n1-1)*(2*pi/N)/pi;%频带宽度figure(2);subplot(2,1,1);plot(x);%录制信号的时域波形title(原始信号的时域波形);%加标题ylabel(幅值/A);%显示纵坐标的表示意义grid;%加网格y0=fft(x);%快速傅立叶变换figure(2);subplot(2,1,2);plot(f,abs(y0(n1);%原始信号的频谱图title(原始信号的频谱图);%加标题xlabel(频率w/pi);%显示横坐标表示的意义ylabel(幅值 );%显示纵坐标表示的意义title(原始信号的频谱图);%加标题grid;%加网格图1.1 原始信号的时域与频谱图1.2滤除无效点针对实际发出声音落后录制动作半拍的现象，如何拔除对无效点的采样的问题: 出现这种现象的原因主要是录音开始时，人的反应慢了半拍，导致出现了一些无效点，而后而出现的无效的点，主要是已经没有声音的动作，先读取声音出来，将原始语音信号时域波形图画出来，根据己得到的信号，可以在第二次读取声音的后面设定采样点，取好有效点，画出滤除无效点后的语音信号时域波形图，对比可以看出。这样就可以解决这个问题。x=wavread(C:UsersacerDesktop数字信号sound.wav, 4000,24000 );%从4000点截取到24000结束plot (x);%画出截取后的时域图形title(截取后的声音时域图形);%标题xlabel(频率);ylabel(振幅) ;grid; %画网格图1.2 去除无效点2.巴特沃斯滤波器的设计2.1设计巴特沃思低通滤波器MATLAB程序如下。滤波器图如图3.3所示。%低通滤波fp=1000;fs=1200;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=buttord(wap,was,rp,rs,s);B,A=butter(N,wc,s);Bz,Az=bilinear(B,A,Fs1);figure(1);h,w=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h);title(巴特沃斯低通滤波器);xlabel(频率（HZ）);ylabel(耗损（dB）);gridon;图2.1 巴特沃思低通滤波器2.2设计巴特沃思高通滤波器MATLAB程序如下。滤波器图如图3.5所示。%高通滤波fp=4800;fs=5000;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=buttord(wap,was,rp,rs,s);B,A=butter(N,wc,high,s);Bz,Az=bilinear(B,A,Fs1);figure(1);h,w=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h);title(巴特沃斯高通滤波器);xlabel(频率（HZ）);ylabel(耗损（dB）);grid on;图2.2巴特沃思高通滤波器2.3设计巴特沃思带通滤波器MATLAB程序如下。滤波器图如图3.7所示。%带通滤波fp=1200,3000;fs=1000,3200;Fs=8000;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);N,wc=buttord(wap,was,rp,rs,s);B,A=butter(N,wc,s);Bz,Az=bilinear(B,A,Fs1);figure(4);h,w=freqz(Bz,Az,512,Fs1*1000);plot(w,abs(h);title(巴特沃斯带通滤波器);xlabel(频率（HZ）);ylabel(耗损（dB）);grid on;图2.3巴特沃思带通滤波器3将声音信号送入滤波器滤波x=wavread(C:UsersacerDesktop数字信号sound.wav);%播放原始信号wavplay(x,fs); %播放原始信号N=length(x);%返回采样点数df=fs/N;%采样间隔n1=1:N/2;f=(n1-1)*(2*pi/N)/pi;%频带宽度figure(4);subplot(4,2,1);plot(x);%录制信号的时域波形title(原始信号的时域波形);%加标题ylabel(幅值/A);%显示纵坐标的表示意义grid;%加网格y0=fft(x);%快速傅立叶变换subplot(4,2,3);plot(f,abs(y0(n1);%原始信号的频谱图title(原始信号的频谱图);%加标题xlabel(频率w/pi);%显示横坐标表示的意义ylabel(幅值 );%显示纵坐标表示的意义title(原始信号的频谱图);%加标题grid;%加网格31低通滤波器滤波fs=8000;beta=10.056;wc=2*pi*1000/fs;ws=2*pi*1200/fs;width=ws-wc;wn=(ws+wc)/2; n=ceil(12.8*pi /width);h=fir1(n,wn/pi,band,kaiser(n+1,beta);h1,w=freqz(h);ys=fftfilt(h,x);%信号送入滤波器滤波，ys为输出fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换figure(4);subplot(4,2,2);%在四行两列的第二个窗口显示图形plot(ys);%信号的时域波形title(低通滤波后信号的时域波形);%加标题xlabel(频率w/pi);ylabel(幅值/A);%显示标表示的意义grid;%网格subplot(4,2,4); %在四行两列的第四个窗口显示图形plot(f, abs(fftwave(n1);%绘制模值xlabel(频率w/pi);ylabel(幅值/A);%显示标表示的意义title(低通滤波器滤波后信号的频谱图);%标题grid;%加网格wavplay(ys,8000); %播放滤波后信号32高通滤波器滤波fs=8000;beta=10.056;ws=2*5000/fs;wc=2*4800/fs;width=ws-wc;wn=(ws+wc)/2; n=ceil(12.8*pi/width);h=fir1(n,wn/pi, high,kaiser(n+2,beta);h1,w=freqz(h);ys=fftfilt(h,x);%将信号送入高通滤波器滤波subplot(4,2,5);%在四行两列的第五个窗口显示图形plot(ys);%信号的时域波形xlabel(频率w/pi);ylabel(幅值/A);%显示标表示的意义title(高通滤波后信号的时域波形);%标题ylabel(幅值/A);%显示纵坐标的表示意义grid;%网格fftwave=fft(ys); %将滤波后的语音信号进行快速傅立叶变换subplot(4,2,7);%在四行两列的第七个窗口显示图形plot(f,abs(fftwave(n1);%绘制模值axis(0 1 0 50);xlabel(频率w/pi);ylabel(幅值/A);%显示标表示的意义title(高通滤波器滤波后信号的频谱图);%标题grid;%加网格wavplay(ys,8000); %播放滤波后信号33带通滤波器fs=8000;beta=10.056;wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;width=ws1-wc1;wn1=(ws1+wc1)/2; wn2=(ws2+wc2)/2;wn=wn1 wn2;n=ceil(12.8/width*pi);h=fir1(n,wn/pi,band,kaiser(n+1,beta);h1,w=freqz(h);ys1= fftfilt(h,x); %将信号送入高通滤波器滤波figure(4);subplot(4,2,6);%在四行两列的第六个窗口显示图形plot(ys1);%绘制后信号的时域的图形title(带通滤波后信号的时域波形);%加标题xlabel(频率w/pi);ylabel(幅值/A);%显示纵坐标表示的意义grid;%网格fftwave=fft(ys1);%对滤波后的信号进行快速傅立叶变换subplot(4,2,8);%在四行两列的第八个窗口显示图形plot(f, abs(fftwave(n1);%绘制模值axis(0 1 0 50);xlabel(频率w/pi);ylabel(幅值/A);%显示标表示的意义title(带通滤波器滤波后信号的频谱图);%加标题grid;%网格wavplay(ys1,8000); %播放滤波后信号图形如下：分析：三个滤波器滤波后的声音与原来的声音都发生了变化。其中低通的滤波后与原来声音没有很大的变化，其它两个都又明显的变化 4.语音信号的回放sound(xlow,Fs,bits); %在Matlab中，函数sound可以对声音进行回放，其调用格式:sound(xhigh, Fs,bits); %sound (x, Fs, bits);sound(xdaitong, Fs,bits);5.男女语音信号的频谱分析5.1 录制一段异性的声音进行频谱分析Fs=8000; %抽样频率time=3; %录音时间fprintf(按Enter键录音%ds,time); %文字提示pause; %暂停命令fprintf(录音中.);x=wavrecord(time*Fs,Fs,double); %录制语音信号fprintf(录音结束); %文字提示fprintf(按Enter键回放录音); pause; %暂停命令wavplay(x,Fs); %按任意键播放语音信号wavwrite(x,Fs,C:UsersacerDesktop数字信号sound2.wav); %存储语音信号5.2 分析男女声音的频谱x=wavread( C:UsersacerDesktop数字信号sound2.wav );%播放原始信号,解决落后半拍wavplay(x,fs); %播放原始信号N=length(x);%返回采样点数df=fs/N;%采样间隔n1=1:N/2;f=(n1-1)*(2*pi/N)/pi;%频带宽度figure(1);subplot(2,2,1);plot(x);%录制信号的时域波形title(原始女生信号的时域波形);%加标题ylabel(幅值/A);%显示纵坐标的表示意义grid;%加网格y0=fft(x);%快速傅立叶变换subplot(2,2,2);plot(f,abs(y0(n1);%原始信号的频谱图title(原始女生信号的频谱图);%加标题xlabel(频率w/pi);%显示横坐标表示的意义ylabel(幅值 );%显示纵坐标表示的意义grid;%加网格y,fs,bits=wavread( C:UsersacerDesktop数字信号sound.wav );% 对语音信号进行采样wavplay(y,fs); %播放原始信号N=length(y);%返回采样点数df=fs/N;%采样间隔n1=1:N/2;f=(n1-1)*(2*pi/N)/pi;%频带宽度subplot(2,2,3);plot(y);%录制信号的时域波形title(原始男生信号的时域波形);%加标题ylabel(幅值/A);%显示纵坐标的表示意义grid;%加网格y0=fft(y);%快速傅立叶变换subplot(2,2,4);%在四行两列的第三个窗口显示图形plot(f,abs(y0(n1);%原始信号的频谱图title(原始男生信号的频谱图);%加标题xlabel(频率w/pi);%显示横坐标表示的意义ylabel(幅值 );%显示纵坐标表示的意义grid;%加网格5.3男女声音的频谱图图5.3男女声音信号波形与频谱对比分析：就时域图看，男生的时域图中振幅比女生的高，对于频谱图女生的高频成分比较多6.噪声的叠加和滤除6.1录制一段背景噪声Fs=8000; %抽样频率time=3; %录音时间fprintf(按Enter键录音%ds,time); %文字提示pause; %暂停命令fprintf(录音中.);x=wavrecord(time*Fs,Fs,double); %录制语音信号fprintf(录音结束); %文字提示fprintf(按Enter键回放录音); pause; %暂停命令wavplay(x,Fs); %按任意键播放语音信号wavwrite(x,Fs,C:UsersacerDesktop数字信号噪音.wav); %存储语音信号6.2 对噪声进行频谱的分析x1,fs,bits=wavread( C:UsersacerDesktop数字信号噪音.wav );%对语音信号进行采样wavplay(x1,fs); %播放噪声信号N=length(x1);%返回采样点数df=fs/N;%采样间隔n1=1:N/2;f=(n1-1)*(2*pi/N)/pi;%频带宽度figure(5);subplot(3,2,1);plot(x1);%信号的时域波形title(噪声信号的时域波形);grid;ylabel(幅值/A);y0=fft(x1);%快速傅立叶变换subplot(3,2,2);plot(f,abs(y0(n1);%噪声信号的频谱图ylabel(幅值);title(噪声信号的频谱图);6.3原始信号与噪音的叠加fs=8000;x,fs,bits=wavread( C:UsersacerDesktop数字信号sound.wav ); %对录入信号进行采样x1,fs,bits=wavread( C:UsersacerDesktop数字信号噪音.wav ); %对噪声信号进行采样yy=x+x1;%将两个声音叠加6.4叠加信号的频谱分析：wavplay(yy,fs); %播放叠加后信号N=length(yy);%返回采样点数df=fs/N;%采样间隔n1=1:N/2;f=(n1-1)*(2*pi/N)/pi;%频带宽度figure(5);subplot(3,2,3);plot(yy,LineWidth,2);%信号的时域波形title(叠加信号的时域波形);xlabel(时间/t);ylabel(幅值/A);grid;y0=fft(yy);%快速傅立叶变换subplot(3,2,4);plot(f,abs(y0(n1);%叠加信号的频谱图title(叠加信号的频谱图);xlabel(频率w/pi);ylabel(幅值/db);grid;6.5 设计一个合适的滤波器将噪声滤除fs=18000;%采样频率Wp=2*1000/fs;%通带截至频率Ws=2*2000/fs;%阻带截至频率Rp=1;%最大衰减Rs=100;%最小衰减N,Wn=buttord(Wp,Ws,Rp,Rs); %buttord函数（n为阶数，Wn为截至频率）num,den=butter(N,Wn); %butter函数（num为分子系数den为分母系数）h,w=freqz(num,den);%DTFT变换ys=filter(num,den,yy);%信号送入滤波器滤波，ys为输出fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换figure(5);subplot(3,2,5);plot(ys);%信号的时域波形title(低通滤波后信号的时域波形);%加标题ylabel(幅值/A);%显示标表示的意义grid;%网格subplot(3,2,6); plot(f, abs(fftwave(n1);%绘制模值title(低通滤波器滤波后信号的频谱图);%标题xlabel(频率w/pi);ylabel(幅值/A);%显示标表示的意义grid;%加网格wavpla