数字信号处理实验六.docx

《数字信号处理实验六.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验六.docx（18页珍藏版）》请在三一办公上搜索。

1、数字信号处理实验六实验六 一、实验名称 离散时间滤波器设计 二、实验目的： 1、掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。 2、加深理解数字滤波器与连续时间滤波器之间的技术转化。 3、掌握脉冲响应不变法设计IIR数字滤波器的优缺点及使用范围。 4、掌握利用双线性变换法设计IIR数字滤波器的原理及具体方法。 5、深入理解利用双线性变换法设计IIR数字滤波器的优缺点及使用范围。 三、实验原理： 1、脉冲响应不变法变换原理 脉冲响应不变法将模拟滤波器的s平面变换成数字滤波器的z平面，从而将模拟滤波器映射成数字滤波器。 IIR滤波器的系数函数为z的有理分式，即 -1H(z)=bzkk

2、=0Nk=1M-k1-akz-k一般满足MN。 转换思路： 变换H(s)拉普拉斯逆变换ha(t)时域采样ha(nT)=h(n)zH(z) 若模拟滤波器的系统函数H只有单阶极点，且假定分母的阶次大于分子的阶次，表达式： H(z)=TAk skT-1zk=11-eNs平面与z平面之间的映射关系。 z=rejws=s+jWz=eresTjw=eesTjWTr=esTw=WT IIR数字滤波器设计的重要环节式模拟低通滤波器的设计，典型的模拟低通滤波器有巴特沃兹和切比雪夫等滤波器。由模拟低通滤波器经过相应的复频率转换为H，由H经过脉冲响应不变法就得到所需要的IIR数字滤波器H。 Matlab信号处理工具

3、箱中提供了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为归一化通带

4、和阻带截止频率; 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

5、,s) 其中：Wp、Ws和Wc均为模拟频率。 切比雪夫I型滤波器设计 切比雪夫I型滤波器为通带波纹控制器：在通带呈现纹波特性，在阻带单调衰减。 N,Wc=cheblord(Wp,Ws,Ap,As) b,a=cheby1(N,Wc,options) 其中参数含义和巴特沃兹的相同。 切比雪夫II型滤波器设计 切比雪夫II型滤波器为阻带波纹控制器：在阻带呈现纹波特性，在通带单调衰减。 N,Wc=cheb2ord(Wp,Ws,Ap,As) b,a=cheby2(N,Wc,options) 其中参数含义和巴特沃兹的相同。 已知模拟滤波器，可以利用脉冲响应不变法转换函数impinvar将其变换为数字滤波器

6、，调用格式为 bz,az=impinvar(b,a,Fs) 其中b，a分别为模拟滤波器系统函数分子、分母多项式系数；Fs为采样频率；bz、az为数字滤波器系统函数的分子、分母多项式系数。设计时要注意模拟原型低通频率预计，否则衰减不能满足设计要求。 5、双线性变换法变换原理 为克服脉冲响应不变法产生频率响应的混叠失真，可以采用非线性频率压缩方法，使s平面与z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，这就是双线性变换法。 转换思路： H(s)写出微分方程近似差分方程写出H(z) 由于双线性变换中，s到z之间的变换时简单的代数关系，得到数字滤波器的系统函数和频率响应，即

7、 1-z-1wwH(z)=Ha(s)|1-z-1=Ha(c) H(e)=H(jW)|=H(jctan) aw-1W=ctans=c1+z221+z-1设模拟系统函数的表达式为 NHa(s)=AskkBk=0k=0NkskA0+A1s+A2s+.+ANsN =B0+B1s+B2s2+.+BNsN应用双线性变换得到H的表达式 H(z)=Ha(s)|s=c1-z-11+z-1=azkN-kbzkk=0k=0N-ka0+a1z-1+a2z-2+.+aNz-N =1+b1z-1+b2z-2+.bNz-Ns平面与z平面之间的映射关系 c+sc+s+jW取模(c+s)2+W2jwz=re=r= 22c-sc

8、-s-jW(c-s)+W用不同的方法选择c可使模拟滤波器频率特性与数字滤波器频率特性在不同频率处有对应的关系。 i、采用使模拟滤波器与数字滤波器在低频处有较确切的对应关系，即在低频处有WW1。当W1较小时，c=2/T. ii、采用数字滤波器的某一特定频率与模拟原型滤波器的一个特定频率Wc严格相对应，则有c=Wccotwc2。 已知模拟滤波器，可以利用双线性变换函数bilinear将其变换为数字滤波器，调用格式为bz,az=bilinear(b,a,Fs) 其中b，a分别为模拟滤波器系统函数分子、分母多项式系数；Fs为采样频率；bz，az为数字滤波器系统函数的分子、分母多项式系数。设计是要注意模

9、拟原型低通频率预畸，否则衰减指标不能满足设计要求。 例3.3用双线性变换设计一个巴特沃斯数字低通滤波器。技术指标为通带截止频率fp=4kHz，通带最大衰减ap=1dB，阻带截止fs=5kHz，阻带最小衰减as=15dB，采样频率fc=30kHz，要求图是滤波器的振幅特性，检验wp,ws对应的衰减指标。 解： clear wp=8*pi*103;ws=10*pi*103;ap=1,as=15; fs=30*103; wp1=wp/fs;ws1=ws/fs; omp1=2*fs*tan(wp1/2);omps=2*fs*tan(ws1/2); N,WC=buttord(omp1,omps,ap,a

10、s,s); b,a=butter(N,WC,s); bz,az=bilinear(b,a,fs); w0=wp1,ws1; Hx=freqz(bz,az,w0); H,W=freqz(bz,az); dbHx=-20*log10(abs(Hx)/max(abs(H); plot(W,abs(H); xlabel(相对频率);ylabel(幅频); grid 1.41.210.80.60.40.2000.511.52相对频率2.533.5四、实验内容 1、要求通带截止频率fp=3kHz，通带最大衰减ap=1dB，阻带截止fs=4.5kHz，阻带最小衰减as=15dB，采样频率fc=30kHz，用

11、脉冲响应不变法设计一个切比雪夫数字低通滤波器，并图示滤波器的振幅特性，检验wp,ws对应的衰减。 2、用双线性变换法设计一个切比雪夫数字滤波器。技术指标为：通带截止频率wp=0.2p，通带最大衰减ap1dB；阻带边缘频率ws=0.3p，阻带最小衰减as15dB。 3.6.2窗函数法设计FIR数字滤波器 1、掌握窗函数法设计FIR数字滤波器的原理即具体方法； 2、深入理解吉布斯现象；理解 不同窗函数的特点。 一、实验原理 1、设计原理 FIR滤波器的设计问题，就是要是所设计的FIR滤波器的频率响应H(e)逼近所要求的理想滤波器的频率响应Hd(e)。逼近可在时域进行，也可以在频域进行。窗函数法设计

12、FIR数字滤波器是在时域进行的，用窗函数截取无限长的hd(n)，这样得到的频率响应jwjwH(ejw)逼近于理想的频率响应Hd(ejw)。 2、设计流程 Hd(ejw)序列傅里叶反变换hd(n)移序加窗截断h(n)序列傅里叶变换H(ejw) 给定希望逼近的频率响应函数Hd(e)； jw1求单位脉冲响应hd(n):(hd(n)=2ppH-pd(ejw)ejwndw) 由过渡带宽及阻带最小衰减的要求，可选定窗形状，并估计窗口长度N。设待求滤波器的过度带用Dw表示，它近似等于窗函数主瓣宽度。因过渡带Dw近似与窗口长度成反比，NA/Dw，A决定与窗口形式； 计算所设计的FIR滤波器的单位脉冲响应； h

13、(n)=hd(n)w(n) 0nN-1 由h(n)求FIR滤波器的频率响应H(e)，检验是否满足设计要求。 一旦选取了窗函数，其指标就是给定的。所以有窗函数设计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,b) %给定beta值的N点凯泽窗函数 jwMatlab中提供的

14、fir1可以用来设计FIR滤波器，调用格式为h=fir1(M,Wc,ftype,window) 其中：h为FIR数字滤波器的系数构成矩阵，Wc是滤波器的截止频率，可以是标量或数组；M+1为FIR数字滤波器的阶数，ftype指定滤波器类型，缺省时为低通，低通用“low”表示，高通用“high”表示，带通用“bandpass”表示，带阻用“stop”表示，window指定窗函数，若不指定，默认为为汉宁窗。 例3.4利用fir1和汉宁窗设计一个N=51，截止频率wc=0.5p的低通滤波器，画出幅频特性。 clear N=51;wc=0.5; h=fir1(50,wc,boxcar(N) H,W=fr

15、eqz(h,1)%数字滤波器频谱数据; plot(W/pi,abs(H); title(矩形窗振幅特性/dB); xlabel(相对频率);ylabel(H(w) 矩形窗振幅特性/dB1.41.210.80.60.40.2000.10.20.30.40.50.6相对频率0.70.80.91例3.5用fir1和布莱克曼窗设计一个N=51，截止频率为wp1=0.3p，wp2=0.4p的带通滤波器。 clear N=51;wc=0.3,0.4; h=fir1(50,wc,bandpass,blackman(N); H,W=freqz(h,1); plot(W/pi,abs(H); title(布莱克

16、曼窗带通振幅特性/dB); xlabel(相对频率);ylabel(H(w) 布莱克曼窗带通振幅特性/dB10.90.80.70.60.50.40.30.20.1000.10.20.30.40.50.6相对频率0.70.80.91二、实验内容 1、窗函数法设计低通数字滤波器， e-jwN/2， 0|w|0.4p Hd(ejw) = 0， 0.4p|w|p N=26，分别利用矩形窗、汉宁窗和布莱克曼窗设计该滤波器，且滤波器具有线性相位。绘出脉冲响应h(n)即滤波器的频率响应； 增加N，观察过渡带和最大肩峰值的变化。 2、利用凯泽窗设计线性相位高通数字滤波器 1， 0.6p|w|p |Hjwd(e

17、)|= 0， 0|w|0.6p 要求N=31，且滤波器具有线性相位。 五、实验结果 3.6.1 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,ws1; Hx=freqz(bz,az,w0); H,W=freqz(bz,az); dbHx=-20*log10(abs(Hx)/max(abs(H) plot(W,abs(H); xlab

18、el();ylabel(); grid bz = 0.0000 0.0054 0.0181 0.0040 0 az = 1.0000 -3.0591 3.8323 -2.2919 0.5495 dbHx = 1.0005 21.5790 图表 1.41.210.80.60.40.2000.511.52相对频率2.533.5结果分析： 2、程序 wp=0.2*pi;ws=0.3*pi;ap=1;as=15; Fs=1; wp1=wp/Fs;ws1=ws/Fs; omp1=2*Fs*tan(wp1/2);oms1=2*Fs*tan(ws1/2); N,WC=cheb2ord(omp1,oms1,a

19、p,as,s) b,a=cheby2(N,as,WC,s); bz,az=bilinear(b,a,Fs) w0=wp1,ws1; Hx=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 = 4 WC = 0.8518 bz = 0.1653 -0.1794 0.2848 -0.1794 0.1653 az = 1.0000 -1.9127 1.7263 -0.6980 0.1408 dbHx = 1.0000 18.2261

20、 图表 10.90.80.70.60.50.40.30.20.1000.511.52相对频率2.533.53.6.2 1、程序 clear N=26;wc=0.4; h1=fir1(25,wc,low,boxcar(N); subplot(321) stem(h1); title(h(n); H1,W1=freqz(h1,1); subplot(322) plot(W1/pi,abs(H1); title(/dB); h2=fir1(25,wc,low,hanning(N); H2,W2=freqz(h2,1); subplot(323) stem(h2); title(h(n); subpl

21、ot(324) plot(W2/pi,abs(H2); title(/dB); h3=fir1(25,wc,low,blackman(N); subplot(325) stem(h3); title(h(n); H3,W3=freqz(h3,1); subplot(326) plot(W3/pi,abs(H3); title(/dB); 图表 矩形窗脉冲响应h(n)0.50-0.50102030汉宁窗脉冲响应h(n)0.50-0.50102030汉宁窗脉冲响应h(n)0.50-0.50102030N=51时 矩形窗低通振幅特性/dB21000.51汉宁窗低通振幅特性/dB21000.51布莱克

22、曼低通振幅特性/dB21000.51当 矩形窗脉冲响应h(n)0.50-0.50204060汉宁窗脉冲响应h(n)0.50-0.50204060汉宁窗脉冲响应h(n)0.50-0.502040602、程序 clear N=31;wc=0.6; h=fir1(30,wc,high,blackman(N); H,W=freqz(h,1); plot(W/pi,abs(H); title(-/dB);xlabel();ylabel(H(w) 图表 矩形窗低通振幅特性/dB21000.51汉宁窗低通振幅特性/dB21000.51布莱克曼低通振幅特性/dB21000.51 凯泽窗高通滤波器振幅特性/dB1.41.210.80.60.40.2000.10.20.30.40.50.6相对频率0.70.80.91六、结果分析 3.6.1 1、脉冲响应不变法dbHx中的1.0005和21.5790为wp,ws处的衰减，可见ap=1.0005(基本 满足)，as=21.5790满足要求 2、双线性变换法dbHx中的1.0000和18.2261为wp,ws处的衰减，可见ap=1.0000，as=18.6621满足要求 3.6.2 观察不同N对应滤波器频率响应波形知，当N增大时，过渡带宽度减小，最大峰值变得 更平稳