《语音信号处理课程设计FIR和IIR数字滤波器设计.doc》由会员分享，可在线阅读，更多相关《语音信号处理课程设计FIR和IIR数字滤波器设计.doc（31页珍藏版）》请在三一办公上搜索。
1、1 绪论数字信号处理是利用计算机或专用处理设备，以数值计算的方法对信号进行采集、抽样、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。它在语音、雷达、图像、系统控制、通信、航空航天、生物医学等众多领域都获得了极其广泛的应用。具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等优点。数字滤波器, 是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲
2、激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。FIR滤波器结构上主要是非递归结构，没有输出到输入的反馈，系统函数H (z)在处收敛，极点全部在z = 0处（因果系统），因而只能用较高的阶数达到高的选择性。FIR数字滤波器的幅频特性精度较之于IIR数字滤波器低，但是线性相位，就是不同频率分量的信号经过fir滤波器后他们的时间差不变，这是很好的性质。FIR数字滤波器是有限的单位响应也有利于对数字信号的处理，便于编程，用于计算的时延也小，这对实时的
3、信号处理很重要。 FIR滤波器因具有系统稳定，易实现相位控制，允许设计多通带（或多阻带）滤波器等优点收到人们的青睐。IIR滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。滤波器的设计可以通过软件或设计专用的硬件两种方式
4、来实现。随着MATLAB软件及信号处理工具箱的不断完善，MATLAB很快成为应用学科等领域不可或缺的基础软件。它可以快速有效地实现数字滤波器的设计、分析和仿真，极大地减轻了工作量,有利于滤波器设计的最优化。2 课程设计内容选择一个语音信号作为分析的对象，或录制一段各人自己的语音信号，对其进行频谱分析；利用MATLAB中的随机函数产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行频谱分析；设计FIR和IIR数字滤波器， 并对被噪声污染的语音信号进行滤波， 分析滤波后信号的时域和频域特征，回放语音信号；最后，设计一个信号处理系统界面。在所设计的系统界面上可以选择滤波器的类型。3 课程设计的具
5、体实现3.1语音信号的采集利用PC 机上的声卡和WINDOWS 操作系统可以进行数字信号的采集。将话筒输入计算机的语音输入插口上,启动录音机。按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音,屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以文件名“speech”保存入g : MATLAB work 中。可以看到,文件存储器的后缀默认为. wav ,这是WINDOWS 操作系统规定的声音文件存的标准。3.2语音信号的时频分析利用MATLAB中的“wavread”命令来读入（采集）语音信号，将它赋值给某一向量。再对其进行采样，记住采样频率和采样点数。下面介绍Wavr
6、ead 函数几种调用格式。（1）y=wavread（file）功能说明：读取file所规定的wav文件，返回采样值放在向量y中。（2）y,fs,nbits=wavread(file) 功能说明：采样值放在向量y中，fs表示采样频率（hz），nbits表示采样位数。（3）y=wavread（file，N）功能说明：读取钱N点的采样值放在向量y中。（4）y=wavread（file，N1,N2）功能说明：读取从N1到N2点的采样值放在向量y中。接下来，对语音信号OriSound.wav进行采样。其程序如下： y,fs,nbits=wavered (OriSound); %把语音信号加载入Matla
7、b 仿真软件平台中然后，画出语音信号的时域波形，再对语音信号进行频谱分析。MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：Xk=fft(xn,N)参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。在本次设计中，我们利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。其程序如下：x,Fs,nbits=wavread(nm.wav);player=audioplayer(x,Fs);play(player);pause(
8、10);stop(player);N=length(x);subplot(2,2,3);plot(x);title(原始信号波形);subplot(2,2,2);Hs=spectrum.welch;psd(Hs,x,Fs,Fs);title(功率图);subplot(2,2,1);specgram(x);title(语谱图);subplot(2,2,4);y=fft(x,N);Xk=abs(y);plot(Xk(1:N/2);title(幅度图);程序结果如下图：图1 原始信号波形及频谱3.3语音信号加噪与频谱分析MATLAB中产生高斯白噪声非常方便，有两个产生高斯白噪声的两个函数 。我们可以
9、直接应用两个函数：一个是WGN，另一个是AWGN。WGN用于产生高斯白噪声，AWGN则用于在某一信号中加入高斯白噪声。也可直接用randn函数产生高斯分布序列。在本次课程设计中，我们是利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。Randn函数有两种基本调用格式：Randn(n)和Randn（m,n）,前者产生nn服从标准高斯分布的随机数矩阵，后者产生mn的随机数矩阵。在这里，我们选用Randn（m,n）函数。语音信号添加噪声及其频谱分析的主要程序如下：y,fs,nbits=wavread (nm.wav);N = lengt
10、h (y) ; %求出语音信号的长度Noise=0.01*randn(size(N); %随机函数产生噪声Si=y+Noise; %语音信号加入噪声 sound(Si);subplot(2,1,1);plot(Si);title(加噪语音信号的时域波形);S=fft(Si); %傅里叶变换subplot(2,1,2);plot(abs(S); title(加噪语音信号的频域波形);程序结果如下图：图2 加噪后的波形及频谱分析3.4 设计FIR和IIR数字滤波器IIR滤波器和FIR滤波器的设计方法完全不同。IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤
11、是：先设计过渡模拟滤波器得到系统函数H（s），然后将H（s）按某种方法转换成数字滤波器的系统函数H(z)。FIR滤波器比鞥采用间接法，常用的方法有窗函数法、频率采样发和切比雪夫等波纹逼近法。对于线性相位滤波器，经常采用FIR滤波器。 对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。具体设计步骤如下：（1）确定所需类型数字滤波器的技术指标。（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为=2/T tan(0.5)(3)将相应类型的模拟滤波器技术
12、指标转换成模拟低通滤波器技术指标。（4）设计模拟低通滤波器。（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。下面我们总结一下利用模拟滤波器设计IIR数字低通滤波器的步骤：（1）确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。（3）按照模拟低通滤波器的技术指标设计及
13、过渡模拟低通滤波器。（4）用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。如前所述，IIR滤波器和FIR滤波器的设计方法有很大的区别。下面我们着重介绍用窗函数法设计FIR滤波器的步骤。如下：（1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。（2）构造希望逼近的频率响应函数。（3）计算h(n).。（4）加窗得到设计结果。接下来，我们根据语音信号的特点给出有关滤波器的技术指标：低通滤波器的性能指标：fp=1000Hz，fc=1
14、200Hz，As=100db ,Ap=1dB高通滤波器的性能指标：fp=3500Hz，fc=4000Hz，As=100dB，Ap=1dB；带通滤波器的性能指标：fp1=1200Hz，fp2=3000hZ，fc1=1000Hz，fc2=3200Hz，As=100dB，Ap=1dB在Matlab中，可以利用函数fir1设计FIR滤波器，利用函数butter,cheby1和ellip设计IIR滤波器，利用Matlab中的函数freqz画出各步步器的频率响应。hn=fir1(M，wc，window)，可以指定窗函数向量window。如果缺省window参数，则fir1默认为哈明窗。其中可选的窗函数有R
15、ectangular Barlrtt Hamming Hann Blackman窗，其相应的都有实现函数。MATLAB信号处理工具箱函数buttp buttor butter是巴特沃斯滤波器设计函数，其有5种调用格式，本课程设计中用到的是N,wc=butter(N,wc,Rp,As,s),该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：B,A=cheby1(N,Rp,wpo,ftypr)B,A=cheby1(N,Rp,wpo,
16、ftypr,s) 函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。下面我们将给出FIR和IIR数字滤波器的主要程序。%=IIR低通滤波器=clear all;Fp=50000;Fs=60000;Ft=200000;wp=2*pi*Fp/Ft; %通带边界频率ws=2*pi*Fs/Ft; %阻带截止频率fp=2*Ft*tan(wp/2);fs=2*Fs*tan(wp/2);n11,wn11=buttord(wp,ws,1,50,s); %求低通滤波器的阶数和截止频率b11,a11=butter(n11,wn11,
17、s); %求S域的频率响应的参数num11,den11=bilinear(b11,a11,0.5); %双线性变换实现S域到Z域的变换h,w=freqz(num11,den11); %根据参数求出频率响应plot(w*8000*0.5/pi,abs(h),g);legend(用butter设计);grid on;图3 IIR低通滤波器%=IIR带通=clear all;Fp1=10000;Fp2=25000;Fs1=10000;Fs2=45000;Ft=200000;wp1=tan(pi*Fp1/Ft); %带通到低通滤波器的转换wp2=tan(pi*Fp2/Ft);ws1=tan(pi*Fs
18、1/Ft);ws2=tan(pi*Fs2/Ft);w=wp1*wp2/ws2;bw=wp2-wp1;wp=1;ws=(wp1*wp2-w.2)/(bw*w);n12,wn12=buttord(wp,ws,1,50,s); %求低通滤波器阶数和截止频率b12,a12=butter(n12,wn12,s); %求S域的频率响应参数num2,den2=lp2bp(b12,a12,sqrt(wp1*wp2),bw);%将S域低通参数转为带通num12,den12=bilinear(num2,den2,0.5); %双线性变换实现S域到Z域的转换h,w=freqz(num12,den12); %根据参数
19、求出频率响应plot(w*8000*0.5/pi,abs(h),g);legend(用butter设计);grid on;图4 IIR带通滤波器%=IIR高通=clear all;Fp=50000;Fs=35000;Ft=200000;wp1=tan(pi*Fp/Ft); %高通到低通滤波器参数转换ws1=tan(pi*Fs/Ft);wp=1;ws=wp1*wp/ws1;n13,wn13=cheb1ord(wp,ws,1,50,s); %求模拟的低通滤波器阶数和截止频率b13,a13=cheby1(n13,1,wn13,s); %求S域的频率响应的参数num,den=lp2hp(b13,a13
20、,wn13); %将S域低通参数转为高通的num13,den13=bilinear(num,den,0.5); %利用双线性变换实现S域到Z域转换h,w=freqz(num13,den13);plot(w*21000*0.5/pi,abs(h),g);title(IIR高通滤波器);legend(用ch1by1设计 );grid on;图5 IIR 高通滤波器%=IIR带阻=clear all;Ft=200000;Fp1=15000;Fp2=35000;Fs1=20000;Fs2=30000;Ap=2;As=13;wp1=2*pi*Fp1/Ft;wp2=2*pi*Fp2/Ft;ws1=2*pi
21、*Fs1/Ft;ws2=2*pi*Fs2/Ft;Wp1=tan(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);Bw=Ws2-Ws1;Wp1=Ws1*Ws2/Wp2;W0=(Ws1*Ws2)0.5;Ws=1;Wp=(Bw*Wp1)/(Ws1*Ws2)-Wp12); N,WN=buttord(Wp,Ws,Ap,As,s); %求模拟的低通滤波器阶数和截止频率B,A=butter(N,WN,s); %求S域的频率响应的参数BT,AT=lp2bs(B,A,W0,Bw); %将S域低通参数转为带阻的num,den=bilinear(BT,AT,0
22、.5); %利用双线性变换实现S域到Z域转换H,W=freqz(num,den,512);plot(W/pi,20*log10(abs(H),g);title(巴特沃兹带阻滤波器);grid on;图6 IIR带阻滤波器%=FIR低通=clear all;wpf=800;wsf=860;Ft=8000;wpz=2*pi*wpf/Ft;wsz=2*pi*wsf/Ft;B=wsz-wpz; %计算过渡带宽度N=ceil(8*pi/B); %求阶wc=(wpz+B/2)/pi; %计算理想高通滤波器通带截止频率hn=fir1(N-1,wc); %调用firl计算低通FIR数字滤波器的h(n)M=80
23、00;hk=fft(hn,M);plot(20*log10(abs(hk),g);title(FIR低通滤波器);grid on;图7 FIR低通滤波器%=FIR高通=clear all;wpf=1260;wsf=1200;Ft=8000;wpz=2*pi*wpf/Ft;wsz=2*pi*wsf/Ft;B=wpz-wsz; %计算过渡带宽度N=ceil(8*pi/B); %求阶wc=(wsz+B/2)/pi; %计算理想高通滤波器通带截止频率hn=fir1(N-1,wc,high); %调用firl计算高通FIR数字滤波器的h(n)M=8000;hk=fft(hn,M);plot(20*log
24、10(abs(hk),g);title(FIR高通滤波器);grid on;图8 FIR高通滤波器%=FIR带通=clear all;wpf=860,2500;wsf=800,2560;Ft=8000; %设置滤波器参数wpz=2*pi*wpf/Ft;wsz=2*pi*wsf/Ft; %模拟频率转化为数字频率B=wpz(1)-wsz(1); %计算过滤带的宽度N0=ceil(8*pi/B);N=N0+mod(N0+1,2); %保证N为奇数wc=(wpz+B/2)/pi;hn=fir1(N-1,wc);M=8000; hk=fft(hn,M); %求h(n)的FFTplot(20*log10(
25、abs(hk),g);title(FIR带通);grid on;图9 FIR带通滤波器%=FIR带阻=clear all;wpf=800,2560;wsf=860,2500;Ft=8000; %设置滤波器参数wpz=2*pi*wpf/Ft;wsz=2*pi*wsf/Ft; %模拟频率转化为数字频率B=wsz(1)-wpz(1); %计算过滤带的宽度N=N0+mod(N0+1,2); %保证N为奇数wc=(wpz+B/2)/pi;hn=fir1(N-1,wc,stop);M=8000; hk=fft(hn,M); %求h(n)的FFTplot(20*log10(abs(hk),g);title(
26、FIR);grid on;图10 FIR带阻滤波器 3.5 用滤波器对加噪语音信号进行滤波用自己设计的各滤波器分别对加噪的语音信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。函数fftfilt用的是重叠相加法实现线性卷积的计算。调用格式为：y=fftfilter(h,x,M)。其中，h是系统单位冲击响应向量；x是输入序列向量；y是系统的输出序列向量；M是有用户选择的输入序列的分段长度，缺省时，默认的输入向量的重长度M=512。函数filter的调用格式：yn=filter(B，A.xn)，它是按照直线型结构实现对
27、xn的滤波。其中xn是输入信号向量，yn输出信号向量。3.6 比较滤波前后语音信号的波形及频谱% =双线性变换法=%*低通滤波器*clear all;Fp=50000;Fs=60000;Ft=200000;wp=2*pi*Fp/Ft; %通带边界频率ws=2*pi*Fs/Ft; %阻带截止频率fp=2*Ft*tan(wp/2);fs=2*Fs*tan(wp/2);n11,wn11=buttord(wp,ws,1,50,s); %求低通滤波器的阶数和截止频率b11,a11=butter(n11,wn11,s); %求s域的频率响应参数num11,den11=bilinear(b11,a11,0.
28、5); %双线变换实现s域到z域的变换y,fs,nbits=wavread (nm.wav);N= length (y) ; %求出语音信号的长度noise=0.01*randn(size(N); %随机函数产生噪声m=y+noise; %语音信号加入噪声M=fft(m); %傅里叶变换z11=filter(num11,den11,m);sound(z11);m11=fft(z11); figure(1); %滤波后的信号subplot(2,2,1);plot(abs(M),g);title(滤波前的信号的频谱);grid on;subplot(2,2,2);plot(abs(m11),r);
29、title(滤波后信号的频谱);grid on;subplot(2,2,3);plot(m);title(滤波前信号的波形);grid on;subplot(2,2,4);plot(z11);title(滤波后信号的波形);grid on;图11 双线性法低通滤波%*高通滤波器*clear all;Fp=50000;Fs=35000;Ft=200000;wp1=tan(pi*Fp/Ft); %通带边界频率ws1=tan(pi*Fs/Ft);wp=1;ws=wp1*wp/ws1;n13,wn13=cheb1ord(wp,ws,1,50,s); %求高通滤波器的阶数和截止频率b13,a13=che
30、by1(n13,1,wn13,s); %求s域的频率响应参数 num,den=lp2hp(b13,a13,wn13); %双线变换实现s域到z域的变换num13,den13=bilinear(num,den,0.5); h,w=freqz(num13,den13);y,fs,nbits=wavread (nm.wav);N= length(y); %求出语音信号的长度Noise=0.01*randn(size(N); %随机函数产生噪声m=y+Noise; %语音信号加入噪声M=fft(m); %傅里叶变换z13=filter(num13,den13,m);sound(z13);m13=fft
31、(z13); %滤波后的信号figure(2);subplot(2,2,1);plot(abs(M),g);title(滤波前的信号的频谱);grid on;subplot(2,2,2);plot(abs(m13),r);title(滤波后信号的频谱);grid on;subplot(2,2,3);plot(m);title(滤波前信号的波形);grid on;subplot(2,2,4);plot(z13);title(滤波后信号的波形);grid on;图12 双线性法高通滤波%*带通滤波器*clear all;Fp1=10000;Fp2=25000;Fs1=10000;Fs2=45000
32、;Ft=200000;wp1=tan(pi*Fp1/Ft); %通带边界频率wp2=tan(pi*Fp2/Ft);ws1=tan(pi*Fs1/Ft);ws2=tan(pi*Fs2/Ft);w=wp1*wp2/ws2;bw=wp2-wp1;wp=1;ws=(wp1*wp2-w.2)/(bw*w);n12,wn12=buttord(wp,ws,1,50,s); %求带通滤波器的阶数和截止频率 b12,a12=butter(n12,wn12,s); num2,den2=lp2bp(b12,a12,sqrt(wp1*wp2),bw);num12,den12=bilinear(num2,den2,0.
33、5); h,w=freqz(num12,den12); %双线变换实现s域到z域的变换 y,fs,nbits=wavread (nm.wav);N=length(y); %求出语音信号的长度Noise=0.01*randn(size(N); %随机函数产生噪声m=y+Noise; %语音信号加入噪声M=fft(m); %傅里叶变换z12=filter(num12,den12,m);sound(z12);m12=fft(z12); %滤波后的信号figure(3);subplot(2,2,1);plot(abs(M),g);title(滤波前的信号的频谱);grid on;subplot(2,2
34、,2);plot(abs(m12),r);title(滤波后信号的频谱);grid on;subplot(2,2,3);plot(m);title(滤波前信号的波形);grid on;subplot(2,2,4);plot(z12);title(滤波后信号的波形);grid on;图13双线性法带通滤波%*带阻滤波器*clear all;Ft=200000;Fp1=15000;Fp2=35000;Fs1=20000;Fs2=30000;Ap=2;As=13;wp1=2*pi*Fp1/Ft;wp2=2*pi*Fp2/Ft;ws1=2*pi*Fs1/Ft;ws2=2*pi*Fs2/Ft;Wp1=t
35、an(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);Bw=Ws2-Ws1;Wp1=Ws1*Ws2/Wp2;W0=(Ws1*Ws2)0.5;Ws=1;Wp=(Bw*Wp1)/(Ws1*Ws2)-Wp12); N,WN=buttord(Wp,Ws,Ap,As,s);B,A=butter(N,WN,s);BT,AT=lp2bs(B,A,W0,Bw);num14,den14=bilinear(BT,AT,0.5);H,W=freqz(num14,den14,512);y,fs,nbits=wavread (nm.wav);N=length(y)
36、; %求出语音信号的长度Noise=0.01*randn(size(N); %随机函数产生噪声m=y+Noise; %语音信号加入噪声M=fft(m); %傅里叶变换z14=filter(num14,den14,m);sound(z14);m14=fft(z14); %滤波后的信号figure(4);subplot(2,2,1);plot(abs(M),g);title(滤波前的信号的频谱);grid on;subplot(2,2,2);plot(abs(m14),r);title(滤波后信号的频谱);grid on;subplot(2,2,3);plot(m);title(滤波前信号的波形);grid on;subplot(2,2,4);plot(z14);title(滤波后信号的波形);grid on;图14双线性法带阻滤波%=窗函数法=%*低通滤波器*x, Fs, nbits=wavread(nm.wav);N=length(x); Noise=0.01*randn(size(N); %产生随机噪声x2=x+Noise;sound(x2,Fs);figure(6);subplot(2,2,1);plot(x2);title(滤波前加噪的信号);grid on;Y=fft(x2,Fs); %
