《语音信号的数字滤波处理没数字信号课程设计.doc》由会员分享，可在线阅读，更多相关《语音信号的数字滤波处理没数字信号课程设计.doc（44页珍藏版）》请在三一办公上搜索。
1、 一语音信号的数字滤波处理 目录1 绪论11.1 课题背景及目的.11.2 设计方案及原理.22 设计题第一部分.62.1 卷积的演示.62.2 采样定理的演示.123 设计题第二部分.153.1 设计要求.153.2 采集语音信号分析.164 设计数字滤波器.194.1 间接法设计IIR数字滤波器.194.1 直接法设计FIR数字滤波器.255 滤波分析.316 心得体会.357 参考文献.35附录.361 绪论1.1 设计目的1.了解语音信号的产生、采集能绘制语音信号的频率响应曲线及频谱图 2. 学会用MATLAB对语音信号进行分析和处理 3. 学习设计滤波器的一般方法，通过课程设计树立正
2、确的设计思想，提高分析问题、解决问题的能力4.掌握用滤波器去除语音信号噪声的方法，观察去噪前后的语音信号。5进行设计基本技能的训练，如查阅设计资料和手册、程序的设计、调试等。1.2 设计要求(1) 熟练掌握数字信号处理的基本概念基本理论。 (2) 熟悉离散信号和离散系统的时域特性。 (3) 掌握序列快速傅里叶变换（FFT）方法。(4) 掌握滤波器的几种网络结构。(5) 学会使用软件MATLAB，掌握其程序设计方法。(6) 掌握利用MATLAB对语音信号进行频谱分析的方法。1.2 课程设计方案选择windows自带的“Windows XP 关机.wav”声音作为语音信号，产生噪声信号加到语音信号
3、中，得到被污染的语音信号，经过频谱分析，分别设计巴特沃思滤波器以及Blackman窗滤波器（低通、带通、高通），对被污染的语音信号滤波，滤掉噪音信号得到符合要求的语音信号。图2.1 总体方案设计流程图图2.2 总体方案设计流程图（Blackman窗）1.2.1 课程设计原理设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数H(z)具有指定的频率特性。数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应，（IIR）数字滤波器和有限长单位冲激响应（FIR）数字滤波器。本次数字信号处理课程设计我们小组选用巴特沃思数字滤波器（IIR）以及Blackman窗数
4、字滤波器（FIR）。1.2.2 IIR数字滤波器的结构与设计IIR滤波器系统函数的极点可以在单位圆内的任何位置，实现IIR滤波器的阶次较低，所用的存储单元较少，效率高，又由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。设计IIR数字滤波器的方法主要有基于冲激响应不变法的IIR数字滤波器设计，基于双线性Z变换法的IIR数字滤波器设计，数字高通、带通及带阻IIR滤波器设计，基于MATLAB函数直接设计IIR数字滤波器。1.2.3 基于双线性Z变换法的IIR数字滤波器设计本次数字信号处理课程设计我们小组选用巴特沃思数字滤波器（IIR）以及Blackman窗数字滤波器（FIR）。由
5、于 的频率映射关系是根据 推导的，所以使j轴每隔2/Ts便映射到单位圆上一周，利用冲激响应不变法设计数字滤波器时可能会导致上述的频域混叠现象。为了克服这一问题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：(1)s平面的整个j轴仅映射为z平面单位圆上的一周。(2)若G(s)是稳定的，由G(s)映射得到的H(z)也应该是稳定的。(3)这种映射是可逆的，既能由G(s)得到H(z)，也能由H(z)得到G(s)。(4)如果G(j0)=1那么，双线性Z变换满足以上4个条件的映射关系，其变换公式为：双线性Z变换的基本思路是，首先将整个s平面压缩到一条从/Ts到/Ts的带宽为2/Ts的横带里，然
6、后通过标准的变换关系将横带变换成整个z平面上去，这样就得到s平面与z平面间的一一对应的单值关系。 在MATLAB中双线性Z变换可以通过bilinear函数实现，其调用格式为Bz,Az=bilinear(B,A,Fs)。 其中B,A为模拟滤波器传递函数G(s)的分子分母多项式的系数向量，而Bz,Az为数字滤波器的传递函数H(z)的分子分母多项式的系数向量。1.2.4 数字高通、带通及带阻IIR滤波器设计除了低通数字滤波器之外，实际中还常常需要高通、带通及带阻数字滤波器。这三种数字滤波器的设计步骤如下：(1)将数字滤波器H(z)的技术指标p和s，通过=tan(/2)转变为模拟滤波器G(s)的技术指
7、标p和s，作归一化处理后，得到p=1，s=s/p；(2)化解为模拟原型滤波器G(s)的技术指标；(3)设计模拟原型滤波器G(p)；(4)将G(p)转换为模拟滤波器的转移函数G(s)；(5)将G(s)转换成数字滤波器的转移函数H(z)s=(z1)(z+1)。 所谓原型滤波器是指归一化的低通滤波器。本节主要讨论通过IIR数字滤波器的原型转换设计法和IIR数字滤波器的直接设计方法来设计数字高通、带通及带阻滤波器其转换方法主要有3种一是直接由模拟低通滤波器转换成数字高通、带通或带阻滤波器二是先由模拟低通滤波器转换成模拟高通、带通或带阻滤波器然后再把它转换成相应的数字滤波器三是将模拟低通滤波器先转换成数
8、字低通滤波器再通过变量代换变换成高通、带通或带阻滤波器。1.3 设计Blackman窗数字滤波器的基本思想1.3.1 FIR数字滤波器的结构与设计FIR滤波器通常采用窗函数方法来设计。正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、海明窗(Hamming window)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartl
9、ett window)及凯塞窗(Kaiser window)。1.3.2 窗设计FIR数字滤波器的基本思想窗设计的基本思想是，首先选择一个适当的理想选频滤波器（它总是具有一个非因果，无限持续时间脉冲响应），然后截取（加窗）它的脉冲响应得到线性相位和因果FIR滤波器。因此这种方法的重点是选择一个合适的窗函数和理想滤波器。我我们用表示理想的选频滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽的低通滤波器由下式给定：为了从得到一个FIR滤波器必须同时在两边截取。而要得到一个因果的线性相位滤波器，它的长度为N必须有：这种操作叫做加窗，h(n)可以看作是与窗函数的乘积：其中根据的不
10、同定义，可以得到不同的窗结构。在频域中，因果FIR滤波器响应由和窗响应的周期卷积得到即：2 设计题第一部分2.1设计要求： 可输入任意两个序列x1(n)、x2(n)，指定x1(n)为自己的学号,例如x1(n)=2,0,1,2,5,7,0,5,0,1,1,4。x2(n)的内容和长度自选。例如x2(n)= 1, 2.43, 6.17,12.93,22.17,32.25,40.88, 45.87, 45.87, 40.88, 32.25, 22.17, 12.93, 6.17, 2.43,1.0000。 分别动态演示两个序列进行线性卷积x1(n)x2(n)和圆周卷积x1(n)x2(n)的过程；要求分
11、别动态演示翻转、移位、乘积、求和的过程； 圆周卷积默认使用两个序列中的最大长度，但卷积前可以指定卷积长度N用以进行混叠分析； 改变圆周卷积长度N，根据实验结果分析2类卷积的关系。 在计算机操作系统下选一段声音文件(XP系统在“C:WINDOWSMedia”)，读取文件取10ms的声音数据产生时域序列x1(n),序列内容自定义。利用x2(n)= 1, 2.43, 6.17,12.93,22.17,32.25,40.88, 45.87, 45.87, 40.88, 32.25, 22.17, 12.93, 6.17, 2.43,1。利用FFT实现快速卷积，验证时域卷积定理，并与直接卷积进行效率对比
12、（验证时采用matlab子函数）。2.2 卷积计算的演示2.2.1 线性卷积的演示演示程序：第 43 页n=-15:18; M=17;yn=zeros(1,19);x1=2,0,1,2,5,7,0,5,0,1,1,4;xm1=zeros(1,11),2,0,1,2,5,7,0,5,0,1,1,4,zeros(1,11);subplot(3,2,1);stem(n,xm1)xlabel(m)ylabel(x1(m)xn2=1, 2.43, 6.17,12.93,22.17,32.25,40.88, 45.87, 45.87, 40.88, 32.25, 22.17, 12.93, 6.17, 2
13、.43,1.0000;subplot(3,2,2);stem(xn2)xlabel(m)ylabel(x2(m)xm2=fliplr(xn2) zeros(1,18);subplot(3,2,3);stem(n,xm2)xlabel(m)ylabel(x2(N-m)title(n=0)yn(1)=sum(xm1.*xm2);for N=1:17xm3=zeros(1,N) fliplr(xn2) zeros(1,M);M=M-1;yn(N+1)=sum(xm1.*xm3);endxm3=zeros(1,18) fliplr(xn2)yn(19)=sum(xm1.*xm3);subplot(3,
14、2,4);stem(yn)xlabel(n)ylabel(y(n)title(分步求卷积)yn1=conv(x1,xn2);subplot(3,2,5);stem(yn1)xlabel(n)ylabel(y(n)title(调用函数直接求)演示结果：2.2.2圆周卷积程序如下：clf;x=2,0,1,2,5,7,0,5,0,1,1,4;y=5,4,3,6,1,7,4,8,4,5,4,8,3,8,7,8;subplot(2,1,1)stem(x,zeros(1,28);N1=length(x);N2=length(y);N=N1;if N2N1 N=N2;endM=input(请输入圆周卷积长度
15、n);if M=N N=M;endif N=N1 x=x;else x=x,zeros(1,N-N1);endX=zeros(N,N);x1=zeros(1,N);pause(2)for n=1:N if n=1 x1(n)=x(n); else x1(n)=x(N-n+2); endendX(1,1:N)=x1;for i=2:N pause(1) for n=1:N if n=1 X(i,n)=X(i-1,N); else X(i,n)=X(i-1,n-1); end end subplot(2,1,1) stem(X(i,1:N);end if N2=N+M-1时，线性卷积与循环卷积相等
16、，且无时域混叠现象。3设计题第二部分（数字滤波器设计）3.1设计要求： 利用Windows下的录音机或其他软件，进行语音信号的采集（*.wav）； 语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图； 产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号； 污染信号的频谱分析，画出被污染的语音信号时域波形和频谱； 根据有关的频谱特性，采用间接法设计IIR数字滤波器，并画出相应滤波器的幅频、相频图（设计3个IIR滤波器）a. 模拟滤波器类型：切比雪夫滤波器（低通、带通、高通）b总体要求：Matlab原程序+仿真波形+技术指标 根据有关的频谱特性，采用直接法设计FIR数字滤波器
17、，并画出相应滤波器的幅频、相频图（设计3个FIR滤波器）a. 滤波器类型：Blackman窗（低通、带通、高通）b总体要求：Matlab原程序+仿真波形+技术指标+窗函数 用自己设计的这些滤波器分别对被不同噪声污染的信号进行滤波； 分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化； 回放语音信号。3.2.1 采集语音信号，频谱分析将电脑D:mediasunshine.wav这一声音文件导入MATLAB软件中，并使用MATLAB绘出采样后的语音信号的时域波形和频谱图。MATLAB程序如下：clf;x,FS,bits=wavread(d:mediasu
18、nshine.wav);x=x(:,1); figure(1); subplot(2,1,1); plot(x); %sound(x,FS,bits); %回放语音 title(语音信号时域波形图)y=fft(x,3260); f=(FS/1630)*1:1630; subplot(2,1,2);plot(f(1:1630),abs(y(1:1630);title(语音信号频谱图); 运行程序得到语音信号时域和频谱图形如图所示:3.2.2添加噪声，频谱分析：MATLAB程序如下：t=0:length(x)-1; zs=0.1*cos(2*pi*10000*t/22050); zs0=0.1*c
19、os(2*pi*10000*t/22050000); figure(2); subplot(2,1,1); plot(zs0) title(噪声信号波形);zs1=fft(zs,1200); subplot(2,1,2);plot(f(1:600),abs(zs1(1:600);title(噪声信号频谱);运行程序得到语音信号时域和频谱图形如图所示:3.2.3添加噪声信号：MATLAB程序如下：y=x+zs;figure(3); subplot(2,1,1); plot(y) title(被噪声污染的信号波形);y1=fft(y,2000); subplot(2,1,2);plot(f(1:1
20、000),abs(y1(1:1000);title(被噪声污染的信号频谱);运行程序得到语音信号时域和频谱图形如图所示:4设计数字滤波器4.1采用间接法设计IIR数字滤波器4.1.1 切比雪夫I型数字低通滤波器先观察信号和噪声的频谱特性，选取合适的通带截止频率和阻带边界频率，故选取wp=2*pi*500;ws=2*pi*2000;4.1.1.1 MATLAB程序如下：wp=2*pi*500;ws=2*pi*1500;Rp=0.1;As=60;N1,wp1=cheb1ord(wp,ws,Rp,As,s);B1,A1=cheby1(N1,Rp,wp1,s);fk=0:1500/512:1500;w
21、k=2*pi*fk;Hk=freqs(B1,A1,wk);plot(fk/1000,20*log(abs(Hk);xlabel(频率（Hz）);ylabel(幅度（dB）)axis(0,1.2,-160,5)4.1.1.2 给被污染的信号加低通数字滤波器wp=2*pi*500;ws=2*pi*2000;Rp=0.1;As=60;N1,wp1=cheb1ord(wp,ws,Rp,As,s);B,A=cheby1(N1,Rp,wp1,s);fk=0:2000/512:2000;wk=2*pi*fk;Hk=freqs(B,A,wk);figure(4);plot(fk/1000,20*log(abs
22、(Hk);xlabel(频率（Hz）);ylabel(幅度（dB）)title(滤波器的频谱分析);axis(0,1.2,-160,5);Bz,Az=impinvar(B,A,Fs); %脉冲响应不变法M=filter(Bz,Az,y);sound(M);%回放污染语音 figure(5);subplot(2,1,1);plot(M);title(滤波后被污染的信号波形);M1=fft(M,2000); subplot(2,1,2);plot(f(1:1000),abs(M1(1:1000);title(滤波后被污染的信号频谱);运行程序得到语音信号时域和频谱图形如图所示:4.1.2设计数字带
23、通滤波器先将带通滤波器的指标装换为对应的低通数字滤波气的指标，设计相应的低通系统函数，对低通系统函数进行频率变换，得到带通数字滤波器4.1.2.1 MATLAB程序：wp=2*pi*100,500;ws=2*pi*10,800;Rp=1;As=20;N1,wp1=cheb1ord(wp,ws,Rp,As,s);B1,A1=cheby1(N1,Rp,wp1,s);fk=0:1500/512:1500;wk=2*pi*fk;Hk=freqs(B1,A1,wk);plot(fk/1000,20*log(abs(Hk);xlabel(频率（Hz）);ylabel(幅度（dB）)axis(0,1.2,-
24、160,5)4.1.2.2 给被污染的信号加带通数字滤波器MATLAB程序：wp=2*pi*100,500;ws=2*pi*10,800;Rp=1;As=20;N1,wp1=cheb1ord(wp,ws,Rp,As,s);B,A=cheby1(N1,Rp,wp1,s);fk=0:2000/512:2000;wk=2*pi*fk;Hk=freqs(B,A,wk);figure(4);plot(fk/1000,20*log(abs(Hk);xlabel(频率（Hz）);ylabel(幅度（dB）)title(带通滤波器的频谱分析);axis(0,1.2,-160,5);Bz,Az=impinvar
25、(B,A,Fs); M=filter(Bz,Az,y);sound(M); %回放污染语音 figure(5);subplot(2,1,1);plot(M);title(滤波后被污染的信号波形);M1=fft(M,2000); subplot(2,1,2);plot(f(1:1000),abs(M1(1:1000);title(滤波后被污染的信号频谱);4.1.2.3运行结果：4.1.3 设计数字高通滤波器4.1.3.1 设计指标：先将高通滤波器的指标装换为对应的低通数字滤波气的指标，设计相应的低通系统函数，对低通系统函数进行频率变换，得到高通数字滤波器，本次取高通滤波器通带边界频率350Hz
26、,阻带边界截止频率150Hz4.1.3.2 MATLAB程序：wp=2*pi*350;ws=2*pi*150;rp=3;rs=10;fs=1000;n,wc=cheb2ord(wp,ws,rp,rs,s);z,p,k=cheb2ap(n,rs);a,b,c,d=zp2ss(z,p,k);at1,bt1,ct1,dt1=lp2hp(a,b,c,d,wc);at2,bt2,ct2,dt2=bilinear(at1,bt1,ct1,dt1,fs);%双线性变换法num,den=ss2tf(at2,bt2,ct2,dt2);figure(4);freqz(num,den,128,fs);axis(15
27、0,400,-60,0)grid on;4.1.3.3 给被污染的信号加高通数字滤波器wp=2*pi*350;ws=2*pi*150;rp=3;rs=10;fs=1000;n,wc=cheb2ord(wp,ws,rp,rs,s);z,p,k=cheb2ap(n,rs);a,b,c,d=zp2ss(z,p,k);at1,bt1,ct1,dt1=lp2hp(a,b,c,d,wc);at2,bt2,ct2,dt2=bilinear(at1,bt1,ct1,dt1,fs);%双线性变换法num,den=ss2tf(at2,bt2,ct2,dt2);figure(4);freqz(num,den,128
28、,fs);axis(150,400,-60,0)grid on;M=filter(num,den,y);sound(M,FS,bits); %回放污染语音 figure(5);subplot(2,1,1);plot(M);title(滤波后被污染的信号波形);M1=fft(M,2000); subplot(2,1,2);plot(f(1:1000),abs(M1(1:1000);axis(0,1400,0,150)title(滤波后被污染的信号频谱);4.1.3.4 程序运行结果：4.2采用直接法设计FIR数字滤波器4.2.1设计Blackman窗-低通：4.2.1.1 设计指标：选取低通滤波
29、器指标，通带边界截止频率500Hz，阻带边界截止频率1000Hz4.2.1.2 Matlab程序：fp=500;fs=1000; wp=2*pi*fp/FS; ws=2*pi*fs/FS; Bt=ws-wp; N0=ceil(11*pi/Bt); N=N0+mod(N0+1,2); wc=(wp+ws)/2/pi; bz,az=fir1(N-1,wc,blackman(N); h,w=freqz(bz,az);figure(4);subplot(2,1,1); plot(bz,az);axis(0,40,-0.2,0.5);title(h(n)波形);subplot(2,1,2); plot(
30、w/pi*FS,20*log(abs(h);axis(0,3500,-400,2);title(低通-损耗函数曲线);4.2.1.3 给被污染的信号加窗：fp=500;fs=1000; wp=2*pi*fp/FS; ws=2*pi*fs/FS; Bt=ws-wp; N0=ceil(11*pi/Bt); N=N0+mod(N0+1,2); wc=(wp+ws)/2/pi; bz,az=fir1(N-1,wc,blackman(N); h,w=freqz(bz,az);figure(4);subplot(2,1,1); plot(bz,az);axis(0,40,-0.2,0.5);title(h
31、(n)波形);subplot(2,1,2); plot(w/pi*FS,20*log(abs(h);axis(0,3500,-400,2);title(低通-损耗函数曲线);M=filter(bz,az,y);sound(M,FS,bits); %播放加窗后语音 figure(5);subplot(2,1,1);plot(M);title(滤波后被污染的信号波形);M1=fft(M,2000); subplot(2,1,2);plot(f(1:500),abs(M1(1:500);title(滤波后被污染的信号频谱);4.2.1.4 程序运行结果：4.2.2设计Blackman窗-带通4.2.
32、2.1 设计指标：选取带通指标参数：通带截止频率100Hz，500Hz，阻带截止频率30Hz,800Hz4.2.2.2 Matlab程序：wlp=100*2*pi/FS;wup=500*2*pi/FS; wls=30*2*pi/FS;wus=1500*2*pi/FS; Bt=wlp-wls; N0=ceil(11*pi/Bt); N=N0+mod(N0+1,2); wc=(wls+wlp)/2/pi,(wus+wup)/2/pi; bz,az=fir1(N-1,wc,blackman(N); h,w=freqz(bz,az);figure(4);subplot(2,1,1); plot(bz,
33、az);axis(100,200,-0.2,0.5);title(h(n)波形);subplot(2,1,2); plot(w/pi*FS,20*log(abs(h);title(带通-损耗函数曲线);4.2.2.3 给被污染的信号加窗：wlp=100*2*pi/FS;wup=500*2*pi/FS; wls=30*2*pi/FS;wus=1500*2*pi/FS; Bt=wlp-wls; N0=ceil(11*pi/Bt); N=N0+mod(N0+1,2); wc=(wls+wlp)/2/pi,(wus+wup)/2/pi; bz,az=fir1(N-1,wc,blackman(N); h
34、,w=freqz(bz,az);figure(4);subplot(2,1,1); plot(bz,az);axis(100,200,-0.2,0.5);title(h(n)波形);subplot(2,1,2); plot(w/pi*FS,20*log(abs(h);title(带通-损耗函数曲线);%axis(0,1,-600,50);M=filter(bz,az,y);sound(M,FS,bits); %播放加窗后语音 figure(5);subplot(2,1,1);plot(M);title(滤波后被污染的信号波形);M1=fft(M,2000); subplot(2,1,2);pl
35、ot(f(1:1000),abs(M1(1:1000);axis(0,2000,0,150);title(滤波后被污染的信号频谱);4.2.2.4 程序运行结果：4.2.3设计Blackman窗-高通：4.2.3.1 设计指标：选取带通指标参数：通带截止频率600Hz阻带截止频率400Hz4.2.3.2 Matlab程序：h,w=freqz(bz,az);figure(4);subplot(2,1,1); plot(bz,az);axis(0,70,-0.2,1);title(h(n)波形);subplot(2,1,2); plot(w/pi*FS,20*log(abs(h);axis(300
36、,1400,-300,100);title(高通-损耗函fp=600;fs=400; wp=2*pi*fp/FS; ws=2*pi*fs/FS; Bt=wp-ws; N0=ceil(11*pi/Bt); N=N0+mod(N0+1,2); wc=(wp+ws)/2/pi; bz,az=fir1(N-1,wc,high,blackman(N); 数曲线);4.2.3.3 给被污染的信号加窗：fp=600;fs=400; wp=2*pi*fp/FS; ws=2*pi*fs/FS; Bt=wp-ws; N0=ceil(11*pi/Bt); N=N0+mod(N0+1,2); wc=(wp+ws)/2/pi; bz,az=fir1(N-1,wc,high,blackman(N); h,w=freqz(bz,az);figure(4);subplot(2,1,1); plot(bz,az);axis(0,70,-0.2,1);title(h(n)波形);subplot(2,1,2);
