语音信号的数字滤波处理(五)

语音信号的数字滤波处理图1所示，程序如下：n=0:50;A=10;a=10*sqrt(2)*pi;w0=10*sqrt(2)*pi;t0=0.001;t1=0.020;t2=0.005;x=A*exp(-a*n*t0).*sin(w0*n*t0);y1=A*exp(-a*n*t1).*sin(w0*n*t1);y2=A*exp(-a*n*t2).*sin(w0*n*t2);subplot(321);stem(n,x,'r'),gridontitle('离散时间信号')xlabel('n');ylabel('幅度')subplot(322);plot(n,x,'g'),gridontitle('连续时间信号')xlabel('t');ylabel('幅度')subplot(323);stem(n,y1),gridontitle('50Hz连续时间信号')xlabel('n');ylabel('幅度')subplot(324);stem(n,y2,'m'),gridontitle('200Hz理想采样信号序列');xlabel('n');ylabel('幅度')subplot(325);stem(n,x,'k'),gridontitle('1000Hz连续时间信号')xlabel('n');ylabel('幅度')t=0:0.001:1;xx=A*exp(-a*n*t0).*sin(w0*n*t0);X=fft(xx);subplot(326)plot(abs(X));xlabel('f/Hz');ylabel('|X(ejw)|')gridon；图2.4采样定理图示3课程设计的基本原理3.1IIR滤波器的设计原理IIR数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法来进行设计，通常采用模拟滤波器原型有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器等。IIR数字滤波器单位取样响应是时宽无限的h(n)，即n0＜n＜∞。数字滤波器按照实现的方法和结构形式分为递归型或非递归型两类。递归型数字滤波器的当前输出y(n)是输入x(n)的当前值和以前各输入值x(n)，x(n–1)，….，及以前各输出值y(n)，y(n–1)，….的函数。一个N阶递归型数字滤波器(IIR滤波器)的差分方程如公式2-1所示：2-1其中，式中的系数至少有一项不为零。说明必须将延时的输出序列进行反馈。递归系统的传统函数定义如公式2-2所示：2-2递归系统的传递函数H(z)在Z平面上不仅有零点，而且有极点。非递归型数字滤波器当前的输出值y(n)仅为当前及以前的输入序列的函数，而与以前的各个输出值无关，因此从结构上看非递归系统没有反馈环路。3.2FIR滤波器的设计原理FIR数字滤波器的单位取样响应是时宽有限的h(n)，即n1≤n＜n2。一个N阶的非递归型数字滤波器(FIR滤波器)的差分方程如公式2-3所示：2-3差分方程式中的系数等于单位取样响应的序列值h(n)，其系统函数H(z)可以表示为公式2-4所示的形式：2-4H(z)是的多项式，因此它的极点只能在Z平面的原点上。FIR滤波器通常采用窗函数方法来设计。窗设计的基本思想是，首先选择一个适当的理想选频滤波器（它总是具有一个非因果，无限持续时间脉冲响应），然后街区（加窗）它的脉冲响应得到线性相位和因果FIR滤波器。我们用表示理想的选频滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽wc<pi的低通滤波器由下式给定：为了从hd(n)得到一个FIR滤波器，必须同时在两边截取hd（n）。而要得到一个因果的线性相位滤波器，它的h(n)长度为N，必须满足公式2-5：2-5h(n)可以看做是hd(n)与窗函数w（n）的乘积，如公式2-6所示：2-6这种操作叫做加窗，其中根据w（n）的不同定义，可以得到不同的窗结构。在频域中，因果FIR滤波器响应由和窗响应的周期卷积得到，如公式2-7所示:2-7常用的窗函数有矩形窗、汉宁（HANNING）窗、海明（矩形）窗、布莱克曼（BLACKMAN）窗、凯泽（KAISER）窗等。4设计程序的调试及运行结果4.1IIR滤波器的设计(1)选择来自于D盘的“ding.wav”声音作为语音信号(2)根据不同的滤波器选择不同的噪声信号加到语音信号中，得到被污染的语音信号。(3)分别设计巴特沃思滤波器（低通、带通、高通）以及Kaiser窗滤波器（低通、带通、高通）对被污染的语音信号滤波，滤掉相应的噪音信号，得到符合要求的语音信号。方案设计流程图如图4.1所示。（图为巴特沃思滤波器，Kaiser窗滤波器与巴特沃思类似）4.1总体方案设计流图4.11巴特沃斯低通滤波器的设计(1)添加高频的噪声信号zs=0.05*cos(2*pi*10000*t/22050)，MATLAB程序如下。噪声信号的时域及幅频图如图3.2所示。%噪声信号zt=0:length(x)-1;zs=0.05*cos(2*pi*10000*t/22050);zs0=0.05*cos(2*pi*10000*t/22050000);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs,1200);%sound(zs,FS,bits);%回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');图4.2噪声信号时域和频谱图(2)设计巴特沃思低通滤波器，MATLAB程序如下。滤波器图如图3.3所示。%低通滤波fp=3000;fs=3500;Fs=22050;rp=1;rs=10;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(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');gridon;图4.3巴特沃思低通滤波器图4.12巴特沃斯高通滤波器的设计(1)添加低频噪声信号zs0=0.05*cos(2*pi*100*t/22050)，MATLAB程序如下。噪声信号的时域及幅频图如图3.4所示。t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs0,1200);%sound(zs,FS,bits);%回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');图4.4噪声信号时域和频谱图(2)设计巴特沃思高通滤波器，MATLAB程序如下。滤波器图如图3.5所示。%高通滤波fp=600;fs=400;Fs=22050;rp=1;rs=10;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(4);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');gridon;图4.5巴特沃思高通滤波器图4.13巴特沃斯高通滤波器的设计(1)添加低频噪声信号zs1=0.05*cos(2*pi*100*t/22050)，MATLAB程序如下。噪声信号的时域及幅频图如图3.6所示。t=0:length(x)-1;zs0=0.05*cos(2*pi*100*t/22050);figure(2);subplot(2,1,1)plot(zs0)title('噪声信号波形');zs1=fft(zs0,1200);%sound(zs,FS,bits);%回放噪音subplot(2,1,2)plot(f(1:600),abs(zs1(1:600)));title('噪声信号频谱');图4.6噪声信号时域和频谱图(2)设计巴特沃思带通滤波器，MATLAB程序如下。滤波器图如图3.7所示。%带通滤波fp=[600,6000];fs=[400,7000];Fs=22050;rp=1;rs=10;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*22050);plot(w,abs(h));title('巴特沃斯带通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');gridon;图4.7巴特沃思带通滤波器图4.2FIR数字滤波器的设计4.21Kaiser窗低通程序的设计技术指标：fr=1000。源程序如下。仿真波形见图4.21。[y,fs,bit]=wavread('1');fr=1000;wr=2*pi*fr;Window=kaiser(8,beta);b=fir1(7,fr/(fs/2),Window);[h,w]=freqz(b,1);subplot(221);plot(w*fs/(2*pi),abs(h));title('低通滤波');y1=filter(b,1,y);y2=fftfilt(b,y);subplot(222);plot(y);title('原始语音信号');subplot(223);plot(y1);title('IIR滤波后语音信号');subplot(224);plot(y2);title('FIR滤波后语音信号');sound(y1,fs,bit);图4.21凯塞窗低通滤波器仿真图4.22Kaiser窗带通滤波器技术指标：fr=1000,fp=3000。源程序如下。仿真波形见图4.22。[y,fs,bit]=wavread('1');fr=1000;fp=3000wr=2*pi*fr;Window=kaiser(16,beta);b=fir1(15,[fr/(fs/2),fp/(fs/2)],Window);[h,w]=freqz(b,1);subplot(221);plot(w*fs/(2*pi),abs(h));title('带通滤波');y1=filter(b,1,y);y2=fftfilt(b,y);subplot(222);plot(y);title('原始语音信号');subplot(223);plot(y1);title('IIR滤波后语音信号');subplot(224);plot(y2);title('FIR滤波后语音信号');sound(y1,fs,bit);图4.22凯塞窗带通滤波器仿真波形4.23Kaiser窗高通滤波器指标：fr=3000。源程序如下。仿真波形见图4.23。[y,fs,bit]=wavread('1');fr=3000;wr=2*pi*fr;Window=kaiser(8,beta);b=fir1(7,fr/(fs/2),'high');[h,w]=freqz(b,1);subplot(221);plot(w*fs/(2*pi),abs(h));title('低通滤波');y1=filter(b,1,y);y2=fftfilt(b,y);subplot(222);plot(y);title('原始语音信号');subplot(223);plot(y1);title('IIR滤波后语音信号');subplot(224);plot(y2);title('FIR滤波后语音信号');sound(y1,fs,bit);图4.23凯塞窗高通滤波器仿真波形4.3回放语音信号经过以上的加噪声处理后，可在Matlab中用函数sound对声音进行回放。其调用格式：sound(y,Fs)，sound(y)和sound(y,Fs,bits)。可以察觉滤波前后的声音有明显的变化。5心得与体会在本次课程设计的过程中，主要遇到的问题是原始音乐信号的选取。有些音乐信号在MATLAB中运行时出错，开始以为自己的操作问题，在跟同学商讨多次更换音乐信号后得到了解决。其次就是在滤波器参数的选取上，开始时不清楚怎么选取，又在课本中找了好久才找到。学习的过程是相互讨论共同进步的，多多讨论课题中遇到的问题，可以巩固我们的知识掌握能力，增加熟练运用度。从本次课程设计的中心来看，课题是希望将数字信号处理技术应用于某一实际