1、课程设计说明书 NO.10语音信号的数字滤波 利用双线性变换法实现IIR数字滤波器的设计一课程设计的目的通过对常用数字滤波器的设计和实现，掌握数字信号处理的工作原理及设计方法；熟悉用双线性变换法设计 IIR 数字滤波器的原理与方法，掌握利用数字滤波器对信号进行滤波的方法，掌握数字滤波器的计算机仿真方法，并能够对设计结果加以分析。二设计方案论证1.IIR数字滤波器设计方法IIR数字滤波器是一种离散时间系统，其系统函数为假设MN，当MN时,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计实际上是求解滤波器的系数和 ，它是数学上的一种逼近问题，即在规定意
2、义上通常采用最小均方误差准那么去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。2.用双线性变换法设计IIR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了
3、频谱混叠现象，映射关系如图1所示。图1双线性变换的映射关系为了将S平面的整个虚轴j压缩到S1平面j1轴上的-/T到/T段上，可以通过以下的正切变换实现 1式中,T仍是采样间隔。当1由-/T经过0变化到/T时，由-经过0变化到+，也即映射了整个j轴。将式1写成将此关系解析延拓到整个S平面和S1平面，令j=s，j1=s1，那么得再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为： (2) (3)式2与式3是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式1与式2的双线性变换符合映射变换应满足的两点要求。首先,把z=ej
4、，可得 (4)即S平面的虚轴映射到Z平面的单位圆。其次，将s=+j代入式4，得因此由此看出，当0时，|z|0时，|z|1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是防止了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个j轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式4所示，重写如下：上式说明，S平面上与Z平面的成非线性的正切关系，如
5、图2所示。由图2看出，在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但当进一步增加时，增长得越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频局部超过折叠频率而混淆到低频局部去的现象，从而消除了频率混叠现象。图2双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式4及图2所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某
6、一常数这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性，不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3所示。图3双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。3.语音信号的采集：（1） 翻开windows系统中的录音机软件，录入自己的声音“数字信号处理。（2） 将音频文件保存“ZJH.wav（3） 翻开音频转换软件，将录制的转换成单声
7、道，采样速率为8KHz。4 语音信号的频谱分析：语音文件“ZJH1.wav频谱分析程序：文件名为zz.mz1=wavread(ZJH1.wav);y1=z1(1:16384);Y1=fft(y1);n=0:16383;plot(n,Y1);gred;图像输出如图4：图4 语音信号频谱分析图5.滤波器的设计：根据以上IIR数字滤波器设计方法及要求，下面运用双线性变换法基于MATLAB设计一个IIR带通滤波器，其中带通的中心频率为p0=1500hz，通带截止频率p1=0 ,p2=3000hz;通带最大衰减p=3dB;阻带最小衰减s=15dB。设计步骤：(1)根据任务,确定性能指标:在设计带通滤波器
8、之前,首先根据工程实际的需要确定滤波器的技术指标:带通滤波器的阻带边界频率关于中心频率p0几何对称，因此ws1=wp0- (ws2-wp0)=300hz通带截止频率wc1=0,wc2=3000hz;阻带截止频率wr1=0，wr2=3300hz；阻带最小衰减s=3dB和通带最大衰减p=15dB;(2)用=2/Txtan(w/2)对带通数字滤波器H(z)的数字边界频率预畸变，得到带通模拟滤波器H(s)的边界频率主要是通带截止频率p1,p2;阻带截止频率s1，s2的转换。Matlab程序：文件名为kk2.mfs=8000;x1=wavread(zjh1.wav);t=0:1/8000 x(size(
9、x1)-1)/8000;Au=0.03;d=Auxcos(2xpix5000xt);x2=x1+d;wp=0.25xpi;ws=0.3xpi;Rp=1;Rs=15;Fs=8000;Ts=1/Fs;wp1=2/Tsxtan(wp/2); %将模拟指标转换成数字指标ws1=2/Tsxtan(ws/2); N,Wn=buttord(wp1,ws1,Rp,Rs,s); %选择滤波器的最小阶数Z,P,K=buttap(N); %创立butterworth模拟滤波器Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用双
10、线性变换法实现模拟滤波器到数字滤波器的转换H,W=freqz(bz,az); %绘制频率响应曲线figure(1)plot(WxFs/(2xpi),abs(H)gridxlabel(频率Hz)ylabel(频率响应幅度)title(Butterworth)f1=filter(bz,az,x2);figure(2)subplot(2,1,1)plot(t,x2) %画出滤波前的时域图title(滤波前的时域波形);subplot(2,1,2)plot(t,f1); %画出滤波后的时域图title(滤波后的时域波形);sound(f1,8000); %播放滤波后的信号wavwrite(f1,fs,
11、k1.wav);F0=fft(f1,1024);f=fsx(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512); %画出滤波前的频谱图title(滤波前的频谱)xlabel(Hz);ylabel(fuzhi);subplot(2,1,2)F1=plot(f,abs(F0(1:512); %画出滤波后的频谱图title(滤波后的频谱)xlabel(Hz);ylabel(fuzhi);滤波器的幅频特性：图5所设计的滤波器幅频特性三设计结果与分析将语音信号进行数字采样并存为z1，对z1叠加15dB的噪音信号
12、存为z2，对信号进行FFT变换，输出语音信号的时域波形图，幅频特性，z1的实际幅频特性。Matlab程序如下：文件名为zjh.ms=wavread(zjh1.wav);fs=8000;L=2floor(log2(length(s);s1=s(1:L);wavwrite(s1,fs,z1.wav);s2=awgn(s1,15);wavwrite(s2,fs,z2.wav);S1=fft(s1);n=0:L-1;subplot(311);plot(n,s1);k=0:L-1;subplot(312);plot(k,abs(S1);deft=fs/L;subplot(313);plot(k(1:L/
13、2)xdeft,abs(S1(1:L/2);输出如图：将原文件截取2s存为z1.wav，将叠加噪音文件存为z2.wav。图6对原始语音进行分析，Matlab程序如下：文件名为kk.mfs=8000;x1=wavread(zjh1.wav); sound(x1,8000); y1=fft(x1,1024); f=fsx(1:512)/1024;figure(1)plot(x1) title(原始语音信号);xlabel(time n);ylabel(fuzhi n);figure(2)freqz(x1) title(频率响应图)figure(3)subplot(2,1,1);plot(abs(y1(1:512) title(原始语音信号FFT频谱)subplot(2,1,2);plot(f,abs(y1(1:512);title(原始语音信号频谱)xlabel(Hz);ylabel(fuzhi);图7原始语音信号图8原始语音信号的频率响应图滤波器的效果分析：信噪比：X=10xlog10(B/N)A是加噪音信号之后通过滤波器，所得信号总能量 。公式为：A=sum(abs(z1).2)A = 45.3256B是没加噪音信号通过滤波器的能量。公式为：B=sum(abs(z2).2)B = 40.1574N是噪音能量。N=A-B=5.1682X是信噪比。
