数字信号处理实验.docx

实验七一、实验目的加深理解IIR数字滤波器的特性，掌握IIR数字滤波器的设计原理与设计方法，以及IIR数字滤波器的应用。二、 实验原理 N阶IIR数字滤波器的系统函数为: IIR数字滤波器的设计主要通过成熟的模拟滤波器设计方法来实现：将数字滤波器设计指标转换为模拟滤波器设计指标，设计出相应的模拟滤波器H(s)，再经过脉冲响应不变法或双线性变换法得到所需的IIR数字滤波器H(z)。 IIR数字滤波器设计的重要环节是模拟原型低通滤波器的设计，主要包括Butterworth、Chebyshev和椭圆等滤波器 MATLAB 信号处理工具箱中提供了IIR滤波器设计的函数。IIR 滤波器阶数选择 buttord -巴特沃斯(Butterworth)滤波器阶数选择。 cheb1ord -切比雪夫(Chebyshev)I 型滤波器阶数选择。 cheb2ord -切比雪夫(Chebyshev)II 型滤波器阶数选择。 ellipord -椭圆(Elliptic)滤波器阶数选择。IIR 滤波器设计 butter -巴特沃斯(Butterworth)滤波器设计 cheby1 -切比雪夫(Chebyshev)I 型滤波器设计 cheby2 -切比雪夫(Chebyshev)II 型滤波器设计 ellip -椭圆(Elliptic)滤波器设计 maxflat -通用的巴特沃斯(Butterworth)低通滤波器设计 yulewalk -Yule-Walker 滤波器设计(直接数字滤波器设计法) 1. Butterworth滤波器设计 Butterworth滤波器是通带、阻带都单调衰减的滤波器。(1)调用buttord函数确定巴特沃斯滤波器的阶数，格式为 N,Wc = buttord(Wp,Ws,Ap,As)输入参数：Ap,As为通带最大衰减和阻带最小衰减，以dB为单位。 Wp,Ws为归一化通带截频和阻带截频，0Wp,Ws1 。输出参数：N为滤波器的阶数；Wc为截频，0 Wc 1。 (2)调用butter函数设计出巴特沃斯滤波器，格式为 b,a = butter(N,Wc,options)输入参数：N和Wc是buttord函数返回的参数，含义见上。 Options=low, high, bandpass, stop, 分别对应低通、高通、带通、带阻，默认情况下为低通或带通。输出参数：b和a为设计出的IIR数字滤波器H(s)的分子多项式和分母多项式的系数矩阵。2. Chebyshev II 型滤波器设计 Chebyshev II 型滤波器为阻带纹波控制器： 在阻带呈现纹波特性。 N,Wc = cheb2ord(Wp, Ws, Ap, As) b,a = cheby2(N,As,Wc,options)3. 椭圆滤波器设计 椭圆滤波器在通阻带都呈现纹波特性。 N,Wc = ellipord(Wp,Ws,Ap,As) b,a = ellip(N,Ap,As,Wc,options)三、实验内容1.信号 ,确定设计指标，实现各种IIR数字滤波器以实现以下信号处理。（1）设计IIR低通滤波器，滤除 的成分。（2）设计IIR高通滤波器，滤除 的成分。（3）设计IIR带通滤波器，滤除 的成分。（4）设计IIR带阻滤波器，滤除 的成分。 要求利用butterord函数求解滤波器的阶数；利用butter函数设计各IIR数字滤波器；画出滤波器的幅度相应和相位响应；给出IIR数字滤波器的系统函数。(1)clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(1/4,1/2,3,60);b,a=butter(N1,Wc,low); freqz(b,a); figure;axis(0,1,-120,0);y=filter(b,a,x);Y=fftshift(fft(y,N);stem(f,abs(Y); (2)clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(2/3,10/24,3,100);b,a=butter(N1,Wc,high); freqz(b,a); figure;axis(0,1,-120,0);y=filter(b,a,x); Y=fftshift(fft(y,N);stem(f,abs(Y);（3）clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(1/8 8/24,1/32 2/3,3,60);b,a=butter(N1,Wc,bandpass); freqz(b,a); figure;axis(0,1,-120,0);y=filter(b,a,x);Y=fftshift(fft(y,N);stem(f,abs(Y); (4)clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(1/32 2/3,0.2 10/27,3,60);b,a=butter(N1,Wc,stop); freqz(b,a); figure;axis(0,1,-120,0);y=filter(b,a,x);Y=fftshift(fft(y,N);stem(f,abs(Y); 2.某带通滤波器的设计指标为 （1）试分别利用巴特沃斯、切比雪夫型、切比雪夫型和椭圆模拟滤波器，通过脉冲响应不变法设计该带通数字滤波器，画出其频率特性，比较设计结果。（2）试分别利用巴特沃斯、切比雪夫型、切比雪夫型和椭圆模拟滤波器，通过双线性变换法设计该带通数字滤波器，画出其频率特性，比较设计结果。（3）分析比较以上设计结果，有何结论？ (1).脉冲响应不变法clear;omegas=0.2 0.72*pi;omegap=0.3 0.6*pi;Ap=1;As=42; Fs=1;ws=omegas*Fs;wp=omegap*Fs;N,Wc=buttord(wp,ws,Ap,As,s);b,a=butter(N,Wc,bandpass,s); % N,Wc=cheb1ord(wp,ws,Ap,As,s);% b,a=cheby1(N,Ap,Wc,bandpass,s);% N,Wc=cheb2ord(wp,ws,Ap,As,s);% b,a=cheby2(N,As,Wc,bandpass,s);% % N,Wc=ellipord(wp,ws,Ap,As,s);% b,a=ellip(N,Ap,As,Wc,bandpass,s);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512); h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm);grid on右图从上到下分别是buttordCheby1 cheby2 ellipord 型滤波器的频率特性(脉冲响应不变法)(2)双线性变换法clear;omegas=0.2 0.72*pi; omegap=0.3 0.6*pi;Ap=1;As=42;Fs=1;T=1/Fs;wp=2*tan(omegap/2)/T;ws=2*tan(omegas/2)/T;% N,Wc=buttord(wp,ws,Ap,As,s);% b,a=butter(N,Wc,bandpass,s);% N,Wc=cheb1ord(wp,ws,Ap,As,s); % b,a=cheby1(N,Ap,Wc,bandpass,s);% % N,Wc=cheb2ord(wp,ws,Ap,As,s);% b,a=cheby2(N,As,Wc,bandpass,s);N,Wc=ellipord(wp,ws,Ap,As,s);b,a=ellip(N,Ap,As,Wc,bandpass,s);numd,dend=bilinear(b,a,Fs);w=linspace(0,pi,512); h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm);grid ontitle(ellipord)（1）利用DTMF的产生原理合成数字“4”，合成公式为： ， 画出其时域波形和频谱图。（2）根据图中的低通和高通滤波器，画出两个滤波器的幅频响应；并利用 对信号 进行滤波；利用fft函数求得输出的频响，绘出幅频响应曲线。（3）设计图中的第一组的四个带通滤波器，画出滤波器的幅频响应；并利用 对图中相应的输入信号进行滤波，在一张图上用不同颜色画出这四条输出曲线。比较其幅值，BP770Hz滤波器的输出最大，因此确定该DTMF的行频为770Hz。(1)clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);stem(k,x);title(时域波形)figure;X=fftshift(fft(x);omega=(-(N-1)/2):(N-1)/2)/N*2*pi;stem(omega,abs(X);title(频谱图) (2)clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=(-(N-1)/2):(N-1)/2)/N*2*pi;%低通wp=2*941/8000;ws=2*1200/8000;Ap=3;As=60;N,wc=buttord(wp,ws,Ap,As);b,a=butter(N,wc);freqz(b,a);figure%滤波输出y=filter(b,a,x);Y=fftshift(fft(y);stem(omega,abs(Y)clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=(-(N-1)/2):(N-1)/2)/N*2*pi;%高通wp=2*1200/8000;ws=2*1000/8000;Ap=3;As=60;N,wc=buttord(wp,ws,Ap,As);b,a=butter(N,wc,high);freqz(b,a);figure%滤波输出y=filter(b,a,x);Y=fftshift(fft(y);stem(omega,abs(Y) (3)clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=(-(N-1)/2):(N-1)/2)/N*2*pi;wp1=692 702./4000;ws1=687 707./4000;wp2=765 775./4000;ws2=760 780./4000;wp3=847 857./4000;ws3=842 862./4000;wp4=936 946./4000;ws4=931 951./4000;Ap=3;As=30;N1,wc1=buttord(wp1,ws1,Ap,As);b1,a1=butter(N1,wc1,bandpass);N2,wc2=buttord(wp2,ws2,Ap,As);b2,a2=butter(N2,wc2,bandpass);N3,wc3=buttord(wp3,ws3,Ap,As);b3,a3=butter(N3,wc3,bandpass);N4,wc4=buttord(wp4,ws4,Ap,As);b4,a4=butter(N4,wc4,bandpass);figure(1)freqz(b1,a1);figure(2)freqz(b2,a2);figure(3)freqz(b3,a3);figure(4)freqz(b4,a4);以上分别是4个带通滤波器的幅频响应.clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=(-(N-1)/2):(N-1)/2)/N*2*pi;wp1=692 702./4000;ws1=687 707./4000;wp2=765 775./4000;ws2=760 780./4000;wp3=847 857./4000;ws3=842 862./4000;wp4=936 946./4000;ws4=931 951./4000;Ap=3;As=20;N1,wc1=buttord(wp1,ws1,Ap,As);b1,a1=butter(N1,wc1,bandpass);N2,wc2=buttord(wp2,ws2,Ap,As);b2,a2=butter(N2,wc2,bandpass);N3,wc3=buttord(wp3,ws3,Ap,As);b3,a3=butter(N3,wc3,bandpass);N4,wc4=buttord(wp4,ws4,Ap,As);b4,a4=butter(N4,wc4,bandpass);y1=filter(b1,a1,x);Y1=abs(fftshift(fft(y1);y2=filter(b2,a2,x);Y2=abs(fftshift(fft(y2);y3=filter(b3,a3,x);Y3=abs(fftshift(fft(y3);y4=filter(b4,a4,x);Y4=abs(fftshift(fft(y4);plot(omega,Y1,-g,omega,Y2,-r,omega,Y3,-b,omega,Y4,-k);由图中可以看出bp770hz滤波器的输出最大，所以该dtmf的行频为770hz（4）clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=(-(N-1)/2):(N-1)/2)/N*2*pi;wp1=1204 1214./4000;ws1=1199 1219./4000;wp2=1331 1341./4000;ws2=1326 1346./4000;wp3=1472 1482./4000;ws3=1467 1487./4000;Ap=3;As=20;N1,wc1=buttord(wp1,ws1,Ap,As);b1,a1=butter(N1,wc1,bandpass);N2,wc2=buttord(wp2,ws2,Ap,As);b2,a2=butter(N2,wc2,bandpass);N3,wc3=buttord(wp3,ws3,Ap,As);b3,a3=butter(N3,wc3,bandpass);y1=filter(b1,a1,x);Y1=abs(fftshift(fft(y1);y2=filter(b2,a2,x);Y2=abs(fftshift(fft(y2);y3=filter(b3,a3,x);Y3=abs(fftshift(fft(y3);plot(omega,Y1,-g,omega,Y2,-r,omega,Y3,-b);图中可以看出bp1209hz滤波器输出最大由(3)(4)可以得出dtmf信号按键为4(5)x=sin(2*pi*770/8000.*k)+sin(2*pi*1477/8000.*k);检验数字6行频输出列频输出得dtmf信号按键为6(6)答:可以,只要取的各个频率的表示颜色不同即可.四、思考题1. 哪些主要因素直接影响IIR数字滤波器的阶数？从工程概念进行定性解释。答：通带截频，阻带截频，通带最大衰减和阻带最大衰减。2. 巴特沃斯、切比雪夫I 型、切比雪夫II 型和椭圆模拟