实验六离散时间滤波器设计分析解析.doc

实验六 离散时间滤波器设计一 实验目的设计和分析一组滤波器，通过得出的结果及函数图象分析离散时间滤波器的一些性能及性质，并获得对设计过程的深入认识。IIR滤波器的设计1掌握利用脉冲响应不变法设计IIR数字滤波器的原理和具体方法；2加深理解数字滤波器与连续时间滤波器之间的技术指标转化；3掌握脉冲响应不变法设计IIR数字滤波器的优缺点及使用范围；4掌握利用双线性变换法设计IIR数字滤波器的原理及具体方法；5深入理解利用双线性变化法设计IIR数字滤波器的优缺点及使用范围；窗函数法设计FIR滤波器1掌握窗函数法设计FIR数字滤波器的原理及具体方法2深入理解吉布斯现象，理解不同窗函数的特点二 实验原理（一）IIR滤波器的设计1、脉冲响应不变法变换原理脉冲响应不变法将模拟滤波器的s平面变换成数字滤波器的z平面，从而将模拟滤波器映射成数字滤波器。IIR数字滤波器的系统函数为（或z）的有理分式，即 一般满足MN1） 转换思路：若模拟滤波器的H(S)只有单阶极点，且假定分母的阶次大于分子的阶次，表达式： 2） s平面与z平面的映射关系 IIR数字滤波器设计的重要环节是低通模拟滤波器的设计，典型的模拟低通滤波器有巴特沃思和切比雪夫滤波器。由模拟低通滤波器经过相应的复频率转换为H（s），由H（s）经过脉冲响应不变法就得到所需要的IIR数字滤波器H（z）。MATLAB信号处理工具箱提供了IIR滤波器设计的函数，常用的函数：IIR滤波器阶数选择Buttord-巴特沃思滤波器阶数选择Cheb1ord-切比雪夫I型滤波器阶数选择Cheb2ord-切比雪夫I型滤波器阶数选择IIR滤波器的设计Butter - 巴特沃思滤波器的设计Cheby1 - 切比雪夫I型滤波器设计Cheby2 - 切比雪夫II型滤波器设计Maxflat - 通用的巴特沃思低通滤波器设计（二）巴特沃思滤波器的设计巴特沃思滤波器是通带，阻带都单调衰减的滤波器（1） 调用buttord函数确定巴特沃思滤波器的阶数，格式 N,Wc=buttord(Wp,Ws,Ap,As)其中：Wp,Ws归一化通带和阻带截止频率（2）调用butter函数设计巴特沃思滤波器，格式 b,a=butter(N,Wc,options)其中：opitions=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为模拟频率（三）切比雪夫I型滤波器的设计切比雪夫I型滤波器为通带波纹控制器：在通带呈现纹波特性，在阻带单调衰减。 N，Wc=cheblord(Wp,Ws,Ap,As) b,a=cheby1(N,As,Wc,options)（四）切比雪夫II型滤波器的设计切比雪夫II型滤波器维族带波纹控控制器：在阻带呈现波纹特性，在通带单调衰减。 N，Wc=che2ord（Wp，Ws，Ap，As） b,a=chebby2(N,As,Wc,option)已知模拟滤波器，可利用脉冲响应不变法转换函数impinvar将其转换为数字滤波器，调用格式为 bz,az=impinvar(b,a,Fs)其中b，a为模拟滤波器的分子分母系数，Fs为采样频率，bz、az为数字滤波器的分子分母多项式系数。（五）双线形变换法变换原理为克服脉冲响应不变法产生频率响应的混叠失真，可以采用非线性频率压缩方法，使S平面与z平面建立一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，这就双线性变换法。已知模拟滤波器，可以利用双线性变换函数bilinear将其变换为数字滤波器，调用格式为bz,az = bilinear(b,a,Fs)其中b，a非别是模拟滤波器系统函数分子，分母多项式系数；Fs为采样频率；bz，az为数字滤波器系统函数的分子分母多项式系数。（六）双线性变换法变换原理为克服脉冲响应不变法产生频率响应的混叠失真，可以采用非线性频率压缩方法，使s平面与z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。1、 转换思路： 由于双线性变换法中，s到z为简单变换关系，得到得到数字滤波器的系统函数和频率响应，即： 设模拟系统函数的表达式为 应用双线性变换得到H(z)的表达式 2、 s与z之间的映射关系 用不同的方法选择c可使模拟滤波器的频率特性与数字滤波器的频率特性在不同频率处有对应关系。（1） 采用是模拟滤波器与数字滤波器在低频处有叫确切的对应关系，即在低频处，当较小时，c=2/T.(2)采用数字滤波器的某一特定频率（例如截止频率）与原型滤波器的一个特定频率相对应，则有。已知模拟滤波器，可以利用双线性变换函数bilinear将其变换为数字滤波器，调用格式为 bz,az=bilinear(b,a,Fs)设计时要注意模拟原型低通频率预畸，否则衰减指标不能满足设计要求。2、 窗函数法设计FIR数字滤波器（1）FIR滤波器的设计问题，就是要使所设计的FIR滤波器的频率响应逼近所要求的理想滤波器的频率响应。逼近可在时域进行，也可在频域进行。窗函数法设计FIR数字滤波器是在时域进行的，用窗函数截取无限长的，这样得到的频率响应逼近于理想的频率响应。（2）设计流程 （1） 给定希望逼近的频率响应函数；（2） 求单位脉冲响应：（3） 由过渡带宽和阻带最小衰减的要求，可选择窗函数，并估计窗口长度N。设待求的滤波器的过渡带用表示，它近似等于窗的主瓣宽度。因过渡带近似与窗口长度成反比，A决定于窗口形式：（4） 计算所设计的FIR滤波器的单位脉冲响应： (5) 由h(n)求FIR滤波器的频率响应，检验是否满足设计要求。 一旦选取了窗函数，其指标（过渡带宽、阻带衰减）就是给定的，所以由窗函数设计FIR数字滤波器就是由阻带衰减指标确定用什么窗，由过渡带宽估计窗函数的长度N。MATLAB提供数种可调用的窗函数，常用的有： Hd=boxcar(N) %N点矩形窗函数 Hd=triang（N） %N点三角窗函数 Hd=hanning（N） %N点汉宁窗函数 Hd=hanming（N） %N点汉明窗函数 Hd=blackman（N） %N点布莱克曼窗函数 Hd=kaiser（N，） %给定值得N点凯泽窗函数Matlab中提供的flir可用来设计FIR滤波器，调用格式：h=flit(M.Wc,ftype,windows)三 实验内容 1、要求通带截止频率fp=3kHz，通带最大衰减=1dB，阻带截止fs=4.5kHz，阻带最小衰减=15dB，采样频率fc=30kHz，用脉冲响应不变法设计一个切比雪夫数字低通滤波器，并图示滤波器的振幅特性，检验wp，ws对应的衰减.设计过程：最终得出 程序：wp=6*pi*103;ws=9*pi*103;ap=1;as=15; %技术指标Fs=30*103; %采样频率wp0=wp/Fs; %数字频率ws0=ws/Fs; %数字频率N,M=cheb1ord(wp,ws,ap,as,s); b,a=cheby1(N,ap,M,s); %调用切比雪夫1型bz,az=impinvar(b,a,Fs) %调用impinvarw0=wp0,ws0; Hx=freqz(bz,az,w0) ; %给Hx赋值H,W=freqz(bz,az); %调用freqz函数dbHx=-20*log10(abs(Hx)/max(abs(H) %赋值plot(W,abs(H); %绘图xlabel(相对频率);ylabel(幅频);grid on %标注 结果： 图1.切比雪夫低通滤波bz = -0.0000 0.0054 0.0181 0.0040 0az = 1.0000 -3.0591 3.8323 -2.2919 0.5495dbHx = 1.0005 21.5790 dbHx中的1.0005及21.5790即在wp、ws处的衰减，可见ap=1.0005略大于1（稍微不符合要求），as=21.579015，满足要求。2、 用双线性变换法设计一个切比雪夫数字1型高通滤波器，技术指标为：采样频率fc=2KHz,通带截止频率fp=700Hz，通带最大衰减ap1dB，阻带边缘频率fs=500Hz，阻带最小衰减as32dB。设计过程; 应用双线性变换得到H(z)的表达式 c=2/T.、预畸。程序：wp=1400*pi;ws=1000*pi;ap=1;as=32; %技术指标Fs=2*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,high,s); %选高通bz,az=bilinear(b,a,Fs)w0=wp1,ws1;Hx=freqz(bz,az,w0); 图2.高通滤波器H,W=freqz(bz,az);dbHx=-20*log10(abs(Hx)/max(abs(H) %求wp1,ws1对应衰减plot(W/pi,abs(H); %绘制幅频特性xlabel(相对频率);ylabel(振频);grid on运行结果：bz = 0.0084 -0.0335 0.0502 -0.0335 0.0084az = 1.0000 2.3741 2.7057 1.5917 0.4103dbHx = 1.0000 33.1098由dbHx =1.0000 33.1098得 ap=11,as=33.109832,符合要求。3、窗函数法设计低通数字滤波器， 设计过程： 由 由求单位脉冲响应： 要求：1、给定矩形窗w（n）=2、 汉宁窗 3、 布莱克曼窗 再由h(n)求FIR滤波器的频率响应，检验是否满足设计要求。（1）N=26，分别用矩形窗、汉宁窗和布莱克曼窗设计该滤波器，且滤波器具有线性相位。绘出脉冲响应h(n)及滤波器的频率响应；程序：N=26;wc=0.4; %确定参数N=26h=fir1(25,wc,boxcar(N) %选取矩形窗H,W=freqz(h,1) subplot(331),plot(h);gridxlabel(n);ylabel(h(n);title(矩形窗时的脉冲响应);subplot(334),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(矩形窗时的振幅特性/dB);subplot(337),plot(W/pi,angle(H),gridxlabel(相对频率);ylabel(相位);title(矩形窗时的相位特性/dB);h=fir1(25,wc,hanning(N) %汉宁窗H,W=freqz(h,1)subplot(332),plot(h);gridxlabel(n);ylabel(h(n);title(汉宁窗时的脉冲响应);subplot(335),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(汉宁窗时的振幅特性/dB);subplot(338),plot(W/pi,angle(H);gridxlabel(相对频率);ylabel(相位);title(矩形窗时的相位特性/dB);h=fir1(25,wc,blackman(N) %blackman窗H,W=freqz(h,1)subplot(333),plot(h);gridxlabel(n);ylabel(h(n);title(布莱克曼窗时的脉冲响应);subplot(336),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(布莱克曼窗时的振幅特性/dB);subplot(339),plot(W/pi,angle(H);gridxlabel(相对频率);ylabel(相位);title(布莱克曼窗时相位特性/dB); 图3.滤波频响（2）增加N，观察过渡带和最大肩峰值的变化。N=70程序：N=70;wc=0.4; %确定参数，选N=70h=fir1(N-1,wc,boxcar(N) %矩形窗H,W=freqz(h,1)subplot(331),plot(h);gridxlabel(n);ylabel(h(n);title(矩形窗时的脉冲响应);subplot(334),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(矩形窗时的振幅特性/dB);subplot(337),plot(W/pi,angle(H),gridxlabel(相对频率);ylabel(相位);title(矩形窗时的相位特性/dB);h=fir1(N-1,wc,hanning(N) %汉宁窗H,W=freqz(h,1)subplot(332),plot(h);gridxlabel(n);ylabel(h(n);title(汉宁窗时的脉冲响应);subplot(335),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(汉宁窗时的振幅特性/dB);subplot(338),plot(W/pi,angle(H);gridxlabel(相对频率);ylabel(相位);title(矩形窗时的相位特性/dB);h=fir1(N-1,wc,blackman(N) %blackman窗H,W=freqz(h,1)subplot(333),plot(h);gridxlabel(n);ylabel(h(n);title(布莱克曼窗时的脉冲响应);subplot(336),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(布莱克曼窗时的振幅特性/dB);subplot(339),plot(W/pi,angle(H);gridxlabel(相对频率);ylabel(相位);title(布莱克曼窗时相位特性/dB);N=100；程序：N=100;wc=0.4; %取N=100，确定参数h=fir1(N-1,wc,boxcar(N) %选矩形窗H,W=freqz(h,1) %数字滤波器频谱数据subplot(331),plot(h);grid % 绘图xlabel(n);ylabel(h(n); %标注title(矩形窗时的脉冲响应);subplot(334),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(矩形窗时的振幅特性/dB);subplot(337),plot(W/pi,angle(H),gridxlabel(相对频率);ylabel(相位);title(矩形窗时的相位特性/dB);h=fir1(N-1,wc,hanning(N) %汉宁窗H,W=freqz(h,1)subplot(332),plot(h);gridxlabel(n);ylabel(h(n);title(汉宁窗时的脉冲响应);subplot(335),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(汉宁窗时的振幅特性/dB);subplot(338),plot(W/pi,angle(H);gridxlabel(相对频率);ylabel(相位);title(矩形窗时的相位特性/dB);h=fir1(N-1,wc,blackman(N) %blackman窗H,W=freqz(h,1)subplot(333),plot(h);gridxlabel(n);ylabel(h(n);title(布莱克曼窗时的脉冲响应);subplot(336),plot(W/pi,abs(H);gridxlabel(相对频率);ylabel(幅度);title(布莱克曼窗时的振幅特性/dB);subplot(339),plot(W/pi,angle(H);gridxlabel(相对频率);yl