数字信号处理实验6离散时间滤波器设计.docx

信号处理实验实验六：离散时间滤波器设计3.6.11、 实验名称：IIR数字滤波器的设计二、实验目的：1、掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。2、 加深理解数字滤波器与连续时间滤波器之间的技术转化。3、 掌握脉冲响应不变法设计IIR数字滤波器的优缺点及使用范围。4、 掌握利用双线性变换法设计IIR数字滤波器的原理及具体方法。5、 深入理解利用双线性变换法设计IIR数字滤波器的优缺点及使用范围。三、实验原理：1、脉冲响应不变法变换原理脉冲响应不变法将模拟滤波器的s平面变换成数字滤波器的z平面，从而将模拟滤波器映射成数字滤波器。IIR滤波器的系数函数为（或z）的有理分式，即一般满足。转换思路：若模拟滤波器的系统函数H（s）只有单阶极点，且假定分母的阶次大于分子的阶次，表达式：s平面与z平面之间的映射关系。 IIR数字滤波器设计的重要环节式模拟低通滤波器的设计，典型的模拟低通滤波器有巴特沃兹和切比雪夫等滤波器。由模拟低通滤波器经过相应的复频率转换为H（s），由H（s）经过脉冲响应不变法就得到所需要的IIR数字滤波器H（z）。Matlab信号处理工具箱中提供了IIR滤波器设计的函数，常用的函数：IIR滤波器阶数选择Buttord-巴特沃兹滤波器阶数选择。Cheb1ord-切比雪夫I型滤波器阶数选择。Cheb2ord-切比雪夫II型滤波器阶数选择。IIR滤波器设计Butter-巴特沃兹滤波器设计。Cheby1-切比雪夫I型滤波器设计。Cheby2-切比雪夫II型滤波器设计。Maxflat-通过的巴特沃兹低通滤波器设计。2、 巴特沃兹滤波器设计巴特沃兹滤波器式通带、阻带都单调衰减的滤波器。调用buttord函数确定巴特沃兹滤波器的阶数，格式N,Wc=buttord(Wp,Ws,Ap,As)其中：Wp，Ws为归一化通带和阻带截止频率;Ap，As为通带最大和最小衰减，单位为dB;N为滤波器阶数，Wc为3dB截止频率，对于带通和带阻滤波器，Wc=W1,W2为矩阵，W1和W2分别为带通的上下截止频率。调用butter函数设计巴特沃兹滤波器，格式b,a=butter(N,Wc,options)其中：options=low，high，bandpass，stop，默认情况下，为低通和带通。b，a为设计出的IIR数字滤波器的分子多项式和分母多项式的系数。注意，利用以上两个函数也可以设计出模拟滤波器，格式为N,Wc=buttord(Wp,Ws,Ap,As,s)b,a=butter(N,Wc,options,s)其中：Wp、Ws和Wc均为模拟频率。3、 切比雪夫I型滤波器设计切比雪夫I型滤波器为通带波纹控制器：在通带呈现纹波特性，在阻带单调衰减。N,Wc=cheblord(Wp,Ws,Ap,As)b,a=cheby1(N,Wc,options)其中参数含义和巴特沃兹的相同。4、 切比雪夫II型滤波器设计切比雪夫II型滤波器为阻带波纹控制器：在阻带呈现纹波特性，在通带单调衰减。N,Wc=cheb2ord(Wp,Ws,Ap,As)b,a=cheby2(N,Wc,options)其中参数含义和巴特沃兹的相同。已知模拟滤波器，可以利用脉冲响应不变法转换函数impinvar将其变换为数字滤波器，调用格式为bz,az=impinvar(b,a,Fs)其中b，a分别为模拟滤波器系统函数分子、分母多项式系数；Fs为采样频率；bz、az为数字滤波器系统函数的分子、分母多项式系数。设计时要注意模拟原型低通频率预计，否则衰减不能满足设计要求。5、 双线性变换法变换原理为克服脉冲响应不变法产生频率响应的混叠失真，可以采用非线性频率压缩方法，使s平面与z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，这就是双线性变换法。转换思路：由于双线性变换中，s到z之间的变换时简单的代数关系，得到数字滤波器的系统函数和频率响应，即 设模拟系统函数的表达式为应用双线性变换得到H（z）的表达式s平面与z平面之间的映射关系用不同的方法选择c可使模拟滤波器频率特性与数字滤波器频率特性在不同频率处有对应的关系。i、采用使模拟滤波器与数字滤波器在低频处有较确切的对应关系，即在低频处有。当较小时，c=2/T.ii、采用数字滤波器的某一特定频率（例如截止频率）与模拟原型滤波器的一个特定频率严格相对应，则有。已知模拟滤波器，可以利用双线性变换函数bilinear将其变换为数字滤波器，调用格式为bz,az=bilinear(b,a,Fs)其中b，a分别为模拟滤波器系统函数分子、分母多项式系数；Fs为采样频率；bz，az为数字滤波器系统函数的分子、分母多项式系数。设计是要注意模拟原型低通频率预畸，否则衰减指标不能满足设计要求。4、 实验内容1、 要求通带截止频率，通带最大衰减，阻带截止，阻带最小衰减，采样频率，用脉冲响应不变法设计一个切比雪夫数字低通滤波器，并图示滤波器的振幅特性，检验对应的衰减。程序：wp=6*pi*103;ws=9*pi*103;ap=1,as=15;Fs=30*103;wp1=wp/Fs;ws1=ws/Fs;N,WC=cheb1ord(wp,ws,ap,as,s);b,a=cheby1(N,ap,WC,s);bz,az=impinvar(b,a,Fs)w0=wp1,ws1Hx=freqz(bz,az,w0)H,W=freqz(bz,az)dbHx=-20*log10(abs(Hx)/max(abs(H)plot(W,abs(H);xlabel(相对频率);ylabel(幅频);grid程序运行结果N = 4WC = 1.8850e+004bz = -0.0000 0.0054 0.0181 0.0040 0az = 1.0000 -3.0591 3.8323 -2.2919 0.5495dbHx =1.0005 21.57902、 用双线性变换法设计一个切比雪夫数字滤波器。技术指标为：通带截止频率，通带最大衰减；阻带边缘频率，阻带最小衰减。程序：wp=0.2*pi*103;ws=0.3*pi*103;ap=1,as=15;Fs=3*103;wp1=wp/Fs;ws1=ws/Fs;omp1=2*Fs*tan(wp1/2);omps=2*Fs*tan(ws1/2);N,WC=cheb1ord(omp1,omps,ap,as,s);b,a=cheby1(N,ap,WC,s);bz,az=bilinear(b,a,Fs)w0=wp1,ws1Hx=freqz(bz,az,w0)H,W=freqz(bz,az)dbHx=-20*log10(abs(Hx)/max(abs(H)plot(W,abs(H);xlabel(相对频率);ylabel(幅频); grid程序运行结果：N = 4WC = 630.6254bz = 1.0e-003 * 0.0268 0.1073 0.1610 0.1073 0.0268az = 1.0000 -3.7583 5.3425 -3.4028 0.8192dbHx = 0.9998 21.7959五、结果分析1、脉冲响应不变法dbHx中的1.0005和21.5790为wp,ws处的衰减，可见ap=1.0005(基本满足)，as=21.5790满足要求2、双线性变换法dbHx中的0.9998和21.7959为wp,ws处的衰减，可见ap=0.9998（基本满足），as=21.7959满足要求3.6.2一、实验名称：窗函数法设计FIR数字滤波器二、实验目的1、 掌握窗函数法设计FIR数字滤波器的原理即具体方法；2、 深入理解吉布斯现象；理解 不同窗函数的特点。三、实验原理1、 设计原理FIR滤波器的设计问题，就是要是所设计的FIR滤波器的频率响应逼近所要求的理想滤波器的频率响应。逼近可在时域进行，也可以在频域进行。窗函数法设计FIR数字滤波器是在时域进行的，用窗函数截取无限长的，这样得到的频率响应逼近于理想的频率响应。2、 设计流程给定希望逼近的频率响应函数；求单位脉冲响应由过渡带宽及阻带最小衰减的要求，可选定窗形状，并估计窗口长度N。设待求滤波器的过度带用表示，它近似等于窗函数主瓣宽度。因过渡带近似与窗口长度成反比，A决定与窗口形式；计算所设计的FIR滤波器的单位脉冲响应； 由h(n)求FIR滤波器的频率响应，检验是否满足设计要求。一旦选取了窗函数，其指标（过度宽度、阻带衰减）就是给定的。所以有窗函数设计FIR滤波器就是由阻带衰减指标确定用什么窗，由过度宽带估计窗函数的长度N。Matlab中提供了数种可以调用的窗函数，常用的有：hd=boxcar(N) %N点矩形窗函数ht=triang(N) %N点三角窗函数hd=hanning(N) %N点汉宁窗hd=hamming(N) %N点汉明窗函数hd=blackman(N) %N点布莱克曼窗hd=kaiser(N,) %给定beta值的N点凯泽窗函数Matlab中提供的fir1可以用来设计FIR滤波器，调用格式为h=fir1(M,Wc,ftype,window)其中：h为FIR数字滤波器的系数构成矩阵（即系统的单位脉冲响应），Wc是滤波器的截止频率（以为单位），可以是标量或数组；M+1为FIR数字滤波器的阶数，ftype指定滤波器类型，缺省时为低通，低通用“low”表示，高通用“high”表示，带通用“bandpass”表示，带阻用“stop”表示，window指定窗函数，若不指定，默认为为汉宁窗。四、实验内容1、 窗函数法设计低通数字滤波器， ， = 0， N=26，分别利用矩形窗、汉宁窗和布莱克曼窗设计该滤波器，且滤波器具有线性相位。绘出脉冲响应h(n)即滤波器的频率响应； 程序：%矩形窗N=26;wc=0.4;nn=0:N-1;h=fir1(N-1,wc,boxcar(N);H,W=freqz(h,1);subplot(311)stem(nn,h);title(矩形窗脉冲响应);xlabel(n);ylabel(h(n);subplot(312)plot(W/pi,abs(H);title(矩形窗振幅特性/dB);xlabel(相对频率);ylabel(H(W) );subplot(313)plot(W/pi,angle(H);title(矩形窗相频特性/dB);xlabel(相对频率);ylabel(H(w);%汉宁窗N=26;wc=0.4;nn=0:N-1;h=fir1(N-1,wc,hanning(N);H,W=freqz(h,1);subplot(311)stem(nn,h);title(汉宁窗脉冲响应);xlabel(n);ylabel(h(n);subplot(312)plot(W/pi,abs(H);title(汉宁窗振幅特性/dB);xlabel(相对频率);ylabel(H(W) );subplot(313)plot(W/pi,angle(H);title(汉宁窗相频特性/dB);xlabel(相对频率);ylabel(H(w);%布莱克曼窗N=26;wc=0.4;nn=0:N-1;h=fir1(N-1,wc,blackman(N);H,W=freqz(h,1);subplot(311)stem(nn,h);title(布莱克曼窗脉冲响应);xlabel(n);ylabel(h(n);subplot(312)plot(W/pi,abs(H);title(布莱克曼窗振幅特性/dB);xlabel(相对频率);ylabel(H(W) );subplot(313)plot(W/pi,angle(H);title(布莱克曼窗相频特性/dB);xlabel(相对频率);ylabel(H(w);图像：增加N，观察过渡带和最大肩峰值的变化。程序：%矩形窗N1=26;wc=0.4;h1=fir1(N1-1,wc,boxcar(N1);H1,W1=freqz(h1,1);subplot(311)plot(W1/pi,abs(H1);title(矩形窗振幅特性/dB N1=26);xlabel(相对频率);ylabel(H1(W) );N2=36;h2=fir1(N2-1,wc,boxcar(N2);H2,W2=freqz(h2,1);subplot(312)plot(W2/pi,abs(H2);title(矩形窗振幅特性/dB N2=36);xlabel(相对频率);ylabel(H2(W) );N3=46;h3=fir1(N3-1,wc,boxcar(N3);H3,W3=freqz(h3,1);subplot(313)plot(W3/pi,abs(H3);title(矩形窗振幅特性/dB N3=46);xlabel(相对频率);ylabel(H3(W) );%汉宁窗N1=26;wc=0.4;h1=fir1(N1-1,wc,hanning(N1);H1,W1=freqz(h1,1);subplot(311)plot(W1/pi,abs(H1);title(汉宁窗振幅特性/dB N1=26);xlabel(相对频率);ylabel(H1(W) );N2=36;h2=fir1(N2-1,wc,hanning(N2);H2,W2=freqz(h2,1);subplot(312)plot(W2/pi,abs(H2);title(汉宁窗振幅特性/dB N2=36);xlabel(相对频率);ylabel(H2(W) );N3=46;h3=fir1(N3-1,wc,hanning(N3);H3,W3=freqz(h3,1);subplot(313)plot(W3/pi,abs(H3);title(汉宁窗振幅特性/dB N3=46);xlabel(相对频率);ylabel(H3(W) );%布莱克曼窗N1=26;wc=0.4;h1=fir1(N1-1,wc,blackman(N1);H1,W1=freqz(h1,1);subplot(311)plot(W1/pi,abs(H1);title(布莱克曼窗振幅特性/dB N1=26);xlabel(相对频率);ylabel(H1(W) );N2=36;h2=fir1(N2