语音信号处理与滤波(数字信号处理课程设计)

《数字信号处理》课程设计报告题目：语音信号处理与滤波专业：通信工程学号：姓名：指导教师：周老师成绩：二〇一二年一月三日一、课程设计目的综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解。二、课程设计基本要求及分组1、对知识点的掌握要求：（1）熟悉离散信号和系统的时域特性。（2）掌握序列快速傅里叶变换FFT方法。（3）学会MATLAB的使用，掌握MATLAB的程序设计方法。（4）利用MATLAB对语音信号进行频谱分析。（5）掌握MATLAB设计FIR和IIR数字滤波器的方法。三、详细设计过程1、问题描述：录制一段自己的语音信号，取不同的数据点对语音信号进行频谱 分析；对所有数据进行插值和抽取处理，改变抽样率再对信号进行频谱分析； 设计FIR和IIR数字滤波器，并对被抽样后的语音信号进行滤波，分析滤波 后信号的时域和频域特征，回放语音信号。2、详细操作步骤和部分运行结果（1）采集语音信号：利用Windows下的录音机，录制一段自己的话音，时间控制在4s左右；然 后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住 采样频率和采样点数。通过使用wavread函数，理解采样频率、采样位数等 概念。wavread函数调用格式如下：y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。[y,fs,nbits]=wavread(file)，采样值放在向量y中，fs表示采样频率(Hz)， bits表示采样位数。y=wavread(file,N)，读取前N点的采样值放在向量y中。y=wavread(file,[N1，N2])，读取从N1点到N2点的采样值放在向量y中。录制并存储语音的程序代码：(函数名'luyin'被主函数‘main’调用)fs=8000;t=4;a=msgbox('请说话...');x=wavrecord(t*fs,fs,1,'Double');%录制语音信号close(a);wavwrite(x,fs,'D:\MATLAB\R2010b\work\数字信号课程设计杜\课程设计.wav');%存储录音信号原始语音时域信号波形图：程序代码：[x1,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计\课程设计.wav');sound(x1,fs,bits);%播放语音figure(1);subplot(1,1,1);plot(x1);%做原始语音信号的时域图形title('原始语音信号作品');xlabel('时间n');ylabel('音量n');运行结果：（2）对语音信号进行频谱分析：在MATLAB中，利用函数fft对信号进行快速傅里叶变换，得到信号的频 谱特性。首先画出语音信号的时域波形，然后对语音信号进行频谱分析。程序代码：[x1,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 忠\课程设计.wav');y1=fft(x1,1024);f=fs*(0:511)/1024;figure(1)y1=fft(x1);%做length(x1)点的FFTy1=fftshift(y1);%平移，是频率中心为0derta_fs=fs/length(x1);%设置频谱的间隔，分辨率subplot(2,1,1);plot([-fs/2:derta_fs:fs/2-derta_fs],abs(y1));%画出原始语音信号的频谱图title('原始语音信号的频谱杜作品');subplot(2,1,2);plot(abs(y1(1:1024)));title('原始语音信号FFT频谱 杜作品');xlabel('点数N');ylabel('幅值原始语音信号频谱和FFT转换后的频谱图（3）设计数字滤波器和画出频率响应：根据语音信号的特点给出有关滤波器的性能指标：低通滤波器性能指标：fp=1000Hz，fc=1200Hz，As=100dB，Ap=1dB；高通滤波器性能指标：fc=4800Hz，fp=5000Hz，As=100dB，Ap=1dB；带通滤波器性能指标：fp1=1200Hz，fp2=3000Hz，fc1=1000Hz，fc2=3200Hz，As=100dB，Ap=1dB。首先用窗函数法（矩形窗(Rectangularwindow)、三角窗(Triangularwindow)、 汉宁窗(Hanningwindow)、海明窗(Hammingwindow)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshevwindow)、巴特里特窗(Bartlettwindow)及凯 塞窗(Kaiserwindow)。）设计上面要求的三种滤波器。在MATLAB中，利用 函数fir1设计FIR滤波器；然后再用双线性变换法设计上面要求的三种滤波 器，利用函数butte、cheby1和ellip设计IIR滤波器；最后，利用MATLAB 中的函数freqz画出各种滤波器的频率响应。（4）IIR滤波器的设计过程IIR低通滤波器：程序代码：()fs=8000;T=1/fs;duration=3;[y,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 杜\课程设计.wav');soundview('D:\MATLAB\R2010b\work\数字信号课程设计杜\课程设 计.wav');N=length(y);t=(1:N)/fs;df=fs/N;n1=1:N/2;f=(n1-1)*df;subplot(2,1,1);plot(y);title('原始信号时域的图形')xlabel('时间/t')ylabel('幅值')subplot(2,1,2);q=fft(y);plot(f,20*log10(abs(q(n1))));title('原始信号频谱图形')xlabel('频率/Hz')ylabel('幅值/db')fs=8000;fp=1000;fst=1200;ap=1;ast=100;wp=2*fp/fs;wst=2*fst/fs;[N,wc]=cheb1ord(wp,wst,ap,ast);[num,den]=cheby1(N,ap,wc,'low');[h,w]=freqz(num,den);subplot(2,1,1);plot(w/pi,abs(h));grid;xlabel('\omega/\pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型低通滤波器的幅频响应');subplot(2,1,2);plot(w/pi,20*log10(abs(h)));grid;xlabel('\omega/\pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型低通滤波器的幅频响应');y0=filter(num,den,y);fftwav=fft(y0);subplot(2,2,1);plot(y);title('原始信号的波形')xlabel('时间/t')ylabel('幅值')subplot(2,2,2);plot(y0);title('通过契比雪夫Ⅰ型低通滤波器后的信号图形')xlabel('时间/t')ylabel('幅值')subplot(2,2,3);q=fft(y);plot(f,20*log10(abs(q(n1))));title('滤波前的信号频谱图')xlabel('频率/Hz')ylabel('幅值/db')subplot(2,2,4);plot(f,20*log10(abs(fftwav(n1))));title('滤波后信号的频谱图')xlabel('频率/Hz')ylabel('幅值/db')wavwrite(y0,fs,'D:\MATLAB\R2010b\work\数字信号课程设计\ 课程设计IIR_LP.wav');[y0,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 杜\课程设计IIR_LP.wav');%语音回放soundview(y0,fs,'课程设计IIR_LP.wav运行结果如下图： IIR高通滤波器设计 程序代码：fs=8000;T=1/fs;duration=3;[y,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 忠\课程设计.wav');soundview('D:\MATLAB\R2010b\work\数字信号课程设计杜艺忠\课程设 计.wav');N=length(y);t=(1:N)/fs;df=fs/N;n1=1:N/2;f=(n1-1)*df;subplot(2,1,1);plot(y);title('原始信号时域的图形')xlabel('时间/t')ylabel('幅值')subplot(2,1,2);q=fft(y);plot(f,20*log10(abs(q(n1))));title('原始信号频谱图形')xlabel('频率/Hz')ylabel('幅值/db')Fs=16000;fp=4000;fst=3800;ap=1;ast=100;wp=2*fp/Fs;wst=2*fst/Fs;[N,wc]=cheb1ord(wp,wst,ap,ast);[num,den]=cheby1(N,ap,wc,'high');[h,w]=freqz(num,den);subplot(2,1,1);plot(w/pi,abs(h));grid;xlabel('\omega/\pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型高通滤波器的幅频响应');subplot(2,1,2);plot(w/pi,20*log10(abs(h)));grid;xlabel('\omega/\pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型高通滤波器的幅频响应');y0=filter(num,den,y);fftwav=fft(y0);subplot(2,2,1);plot(y);title('原始信号的波形')xlabel('时间/t')ylabel('幅值')subplot(2,2,2);plot(y0);title('通过契比雪夫Ⅰ型高通滤波器后的信号图形')xlabel('时间/t')ylabel('幅值')subplot(2,2,3);q=fft(y);plot(f,20*log10(abs(q(n1))));title('滤波前的信号频谱图')xlabel('频率/Hz')ylabel('幅值/db')subplot(2,2,4);plot(f,20*log10(abs(fftwav(n1))));title('滤波后信号的频谱图')xlabel('频率/Hz')ylabel('幅值/db')wavwrite(y0,fs,'D:\MATLAB\R2010b\work\数字信号课程设计\课程 设计IIR_HP.wav');[y0,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 忠\课程设计IIR_HP.wav');%语音回放soundview(y0,fs,'课程设计IIR_HP.wav');程序运行结果如下图：IIR带通滤波器的设计程序代码:fs=8000;T=1/fs;duration=3;[y,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 忠\课程设计.wav');soundview('D:\MATLAB\R2010b\work\数字信号课程设计杜\课程设 计.wav');N=length(y);t=(1:N)/fs;df=fs/N;n1=1:N/2;f=(n1-1)*df;subplot(2,1,1);plot(y);title('原始信号时域的图形')xlabel('时间/t')ylabel('幅值')subplot(2,1,2);q=fft(y);plot(f,20*log10(abs(q(n1))));title('原始信号频谱图形')xlabel('频率/Hz')ylabel('幅值/db')fs=8000;ast=100;ap=1;wp=[12003000]*2/fs;wst=[10003200]*2/fs;[N,wn]=cheb1ord(wp,wst,ap,ast);[num,den]=cheby1(N,ap,wn);[h,w]=freqz(num,den);subplot(2,1,1);plot(w/pi,abs(h));grid;xlabel('\omega/\pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型带通滤波器的幅频响应');subplot(2,1,2);plot(w/pi,20*log10(abs(h)));grid;xlabel('\omega/\pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型带通滤波器的幅频响应');y1=filter(num,den,y);fftwav=fft(y1);subplot(2,2,1);plot(y);title('原始信号的波形')xlabel('时间/t')ylabel('幅值')subplot(2,2,2);plot(y1);title('通过IIR带通滤波器后的信号图形')xlabel('时间/t')ylabel('幅值')subplot(2,2,3);q=fft(y);plot(f,20*log10(abs(q(n1))));title('滤波前的信号频谱图')xlabel('频率/Hz')ylabel('幅值/db')subplot(2,2,4);plot(f,20*log10(abs(fftwav(n1))));title('带通IIR滤波器滤波后信号的频谱图')xlabel('频率/Hz')ylabel('幅值/db')wavwrite(y1,fs,'D:\MATLAB\R2010b\work\数字信号课程设计杜\课程 设计IIR_BP.wav');[y1,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 忠\课程设计IIR_BP.wav');%语音回放)soundview(y1,fs,'课程设计IIR_BP.wav');程序运行结果：（5）FIR滤波器的设计FIR低通滤波器：程序代码：fs=8000;T=1/fs;duration=3;[y,fs,sbits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 杜艺\课程设计.wav');soundview('D:\MATLAB\R2010b\work\数字信号课程设计杜\课程设 计.wav');N=length(y);t=(1:N)/fs;df=fs/N;n1=1:N/2;f=(n1-1)*df;subplot(2,1,1);plot(y);title('原始信号时域的图形')xlabel('时间/t')ylabel('幅值')subplot(2,1,2);q=fft(y);plot(f,20*log10(abs(q(n1))));title('原始信号频谱图形')xlabel('频率/Hz')ylabel('幅值/db')fs=8000;fp=1000;fb=1200;ap=1;ast=100;wp=2*pi*fp/fs;wst=2*pi*fb/fs;deltaW=wst-wp;wc=(wst+wp)/2;N=ceil(2*6.4*pi/deltaW);hdWindow=ideallp(wc,N);%---ideallp（）函数（非系统自有函数）在系 统安装目录的WORK子目录中，ideallp.mwdWindow=Kaiser(N,10.056);hr=wdWindow.*hdWindow';n=0:N-1;subplot(2,2,1);stem(n,wdWindow);xlabel('时间')ylabel('振幅作品')title('凯泽窗')subplot(2,2,2);stem(n,hr);xlabel('时间')ylabel('振幅')title('凯泽窗冲激响应')[H,W]=freqz(hr,1);subplot(2,2,3);plot(W/pi,abs(H));xlabel('\omega/\pi')ylabel('绝对振幅')title('低通滤波器幅频特性（1）')subplot(2,2,4);plot(W/pi,20*log10(abs(H)));xlabel('\omega/\pi')ylabel('相对振幅')title('低通滤波器幅频特性（2）')y0=fftfilt(hr,y);fftwav=fft(y0);subplot(2,2,1);plot(y);title('原始信号的波形')xlabel('时间/t')ylabel('幅值')subplot(2,2,2);plot(y0);title('通过凯泽低通滤波器后的信号图形')xlabel('时间/t')ylabel('幅值')subplot(2,2,3);q=fft(y);plot(f,20*log10(abs(q(n1))));title('滤波前的信号频谱图')xlabel('频率/Hz')ylabel('幅值/db')subplot(2,2,4);plot(f,20*log10(abs(fftwav(n1))));title('滤波后信号的频谱图')xlabel('频率/Hz')ylabel('幅值/db')wavwrite(y0,fs,'D:\MATLAB\R2010b\work\数字信号课程设计艺\ 课程设计FIR_LP.wav');[y0,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 杜\课程设计FIR_LP.wav');%语音回放soundview(y0,fs,'课程设计FIR_LP.wav');程序运行结果：FIR高通滤波器：程序代码：fs=8000;T=1/fs;duration=3;[y,fs,sbits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计忠\课程设计.wav');soundview('D:\MATLAB\R2010b\work\数字信号课程设计杜\课程设计.wav');N=length(y);t=(1:N)/fs;df=fs/N;n1=1:N/2;f=(n1-1)*df;subplot(2,1,1);plot(y);title('原始信号时域的图形')xlabel('时间/t')ylabel('幅值')subplot(2,1,2);q=fft(y);plot(f,20*log10(abs(q(n1))));title('原始信号频谱图形')xlabel('频率/Hz')ylabel('幅值/db')Fs=16000;fp=4000;fst=3800;ast=100;ap=1;wp=2*pi*fp/Fs;wst=2*pi*fst/Fs;wc=(wst+wp)/2;deltaW=wp-wst;N0=ceil(2*6.4*pi/deltaW);%---查表7-3（P342）凯泽窗N=N0+mod(N0+1,2);%---确保N为奇数N=ceil(2*6.4*pi/deltaW);hdWindow=ideallp(pi,N)-ideallp(wc,N);%---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录中，ideallp.mwdWindow=Kaiser(N,10.056);hr=wdWindow.*hdWindow';n=0:N-1;subplot(2,2,1);stem(n,wdWindow);xlabel('时间')ylabel('振幅')title('凯泽窗')subplot(2,2,2);stem(n,hr);xlabel('时间')ylabel('振幅')title('凯泽窗冲激响应')[H,W]=freqz(hr,1);subplot(2,2,3);plot(W/pi,abs(H));xlabel('\omega/\pi')ylabel('绝对振幅')title('凯泽高通滤波器幅频特性（1）')subplot(2,2,4);plot(W/pi,20*log10(abs(H)));xlabel('\omega/\pi')ylabel('相对振幅')title('凯泽高通滤波器幅频特性（2）')y0=fftfilt(hr,y);fftwav=fft(y0);subplot(2,2,1);plot(y);title('原始信号的波形')xlabel('时间/t')ylabel('幅值')subplot(2,2,2);plot(y0);title('通过凯泽高通滤波器后的信号图形')xlabel('时间/t')ylabel('幅值')subplot(2,2,3);q=fft(y);plot(f,20*log10(abs(q(n1))));title('滤波前的信号频谱图')xlabel('频率/Hz')ylabel('幅值/db')subplot(2,2,4);plot(f,20*log10(abs(fftwav(n1))));title('凯泽高通滤波器滤波后信号的频谱图')xlabel('频率/Hz')ylabel('幅值/db')wavwrite(y0,fs,'D:\MATLAB\R2010b\work\数字信号课程设计杜\课程 设计FIR_HP.wav');[y0,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计 忠\课程设计FIR_HP.wav');%语音回放soundview(y0,fs,'课程设计FIR_HP.wav');程序运行结果：FIR带通滤波器：程序代码：fs=8000;T=1/fs;duration=3;[y,fs,sbits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计\课程设计.wav');soundview('D:\MATLAB\R2010b\work\数字信号课程设计杜\课程设计.wav');N=length(y);t=(1:N)/fs;df=fs/N;n1=1:N/2;f=(n1-1)*df;subplot(2,1,1);plot(y);title('原始信号时域的图形')xlabel('时间/t')ylabel('幅值')subplot(2,1,2);q=fft(y);plot(f,20*log10(abs(q(n1))));title('原始信号频谱图形')xlabel('频率/Hz')ylabel('幅值/db')fs=8000;fp1=1200;fp2=3000;fst1=1000;fst2=3200;ast=100;ap=1;wp1=2*pi*fp1/fs;wp2=2*pi*fp2/fs;wst1=2*pi*fst1/fs;wst2=2*pi*fst2/fs;wc1=(wst1+wp1)/2;wc2=(wst2+wp2)/2;deltaW=min((wp1-wst1),(wst2-wp2));N0=ceil(2*6.4*pi/deltaW);%---查表7-3（P342）凯泽窗N=N0+mod(N0+1,2);%---确保N为奇数N=ceil(2*6.4*pi/deltaW);hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%---ideallp（）函数（非系统 自有函数）在系统安装目录的WORK子目录中，ideallp.mwdWindow=Kaiser(N,10.056);hr=wdWindow.*hdWindow';n=0:N-1;subplot(2,2,1);stem(n,wdWindow);xlabel('时间')ylabel('振幅')title('凯泽窗')subplot(2,2,2);stem(n,hr);xlabel('时间')ylabel('振幅')title('凯泽窗冲激响应')[H,W]=freqz(hr,1);subplot(2,2,3);plot(W/pi,abs(H));xlabel('\omega/\pi')ylabel('绝对振幅')title('带通滤波器幅频特性（1）')subplot(2,2,4);plot(W/pi,20*log10(abs(H)));xlabel('\omega/\pi')ylabel('相对振幅')title('带通滤波器幅频特性（2）')y0=fftfilt(hr,y);fftwav=fft(y0);subplot(2,2,1);plot(y);title('原始信号的波形')xlabel('时间/t')ylabel('幅值')subplot(2,2,2);plot(y0);title('通过带通滤波器后的信号图形')xlabel('时间/t')ylabel('幅值')subplot(2,2,3);q=fft(y);plot(f,20*log10(abs(q(n1))));title('滤波前的信号频谱图')xlabel('频率/Hz')ylabel('幅值/db')subplot(2,2,4);plot(f,20*log10(abs(fftwav(n1))));title('带通滤波器滤波后信号的频谱图')xlabel('频率/Hz')ylabel('幅值/db')wavwrite(y0,fs,'D:\MATLAB\R2010b\work\数字信号课程设计\课程 设计FIR_BP.wav');[y0,fs,bits]=wavread('D:\MATLAB\R2010b\work\数字信号课程设计杜艺 忠\课程设计FIR_BP.wav');%语音回放soundview(y0,fs,'课程设计FIR_BP.wav');程序运行结果：（6）设计系统界面在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤 波器的频率响应，选择信号等。主函数程序（main）：b=menu('忠作品版权所有，翻版必究','录音','原始信号采样后时域图和频谱 图','FIR滤波器','IIR滤波器','退出');ifb==5b=0;Endifb==1luyin;main;elseifb==2temp=menu('作品版权所有，翻版必究','播放原始语音并生成时域图',' 原始语音频谱图及做FFT变换后频谱图','返回');iftemp==1cjyy;main; elseiftemp==2