数字信号处理_IIR及FIR设计

1、IIR数字滤波器设计作业通信工程 2013040157.12设计一个数字切比雪夫I型带通滤波器，给定指标为：(1) 波纹，当 (2) 衰减, ， (3) 抽样频率试用冲激响应不变法，双线性变换法进行设计，最后写出的表达式，并画出系统的幅频响应特性。解：冲激响应不变法设计程序如下%冲激响应不变法，ex712.mclc;clear allOmegaP1=2*pi*200;OmegaP2=2*pi*400;%带通截止频率OmegaS1=2*pi*100;OmegaS2=2*pi*600;%1.5kHz阻带Rp=2;%波纹系数As=20;%阻带衰减dBFs=2*103;%抽样频率2khzOmegaP=

2、OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs;ws=OmegaS/Fs;%等效数字频率N,OmegaC=cheb1ord(OmegaP,OmegaS,Rp,As,s)%滤波器阶数截止频率b,a=cheby1(N,Rp,OmegaC,s);%AF系统函数的分子 分母bz,az=impinvar(b,a,Fs)%冲击不变法AF to DFw0=wp,ws%四个频点Hx=freqz(bz,az,w0);%检验H,w=freqz(bz,az);%计算0pi上的响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化并求dB

3、dbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应plot(w/pi,dbH);%画图xlabel(Omega/pi);ylabel(dB);axis(0,1,-60,5);grid运行,得N = 3;OmegaC =1.0e+03 * 1.2566 2.5133w0 = 0.6283 1.2566 0.3142 1.8850bz = -0.0000 0.0272 -0.0581 0.0109 0.0437 -0.0237 0az = 1.0000 -3.3030 6.0060 -6.7463 5.1356 -2.4093 0.6290dbHx = 2.0022

4、2.0015 41.6739 30.7707图 1冲激响应不变法设计IIR带通滤波器由程序返回得到的数值可以得知，这是一个3阶带通滤波器。转换为模拟原型滤波器的频率见OmegaC返回的值。数字滤波器的边界频率见w0返回的值。由bz,az返回值可得由dbHx返回值可知在四个边界频率处幅度响应大小:对应到模拟频率，即有2004001006002.00222.001541.673930.7707满足题目要求的设计指标，通带波纹不大于2dB, 截止频率处衰减不小于20dB。双线性变换法设计程序如下%双线性变换法，ex7122mclc;clear allOmegaP1=2*pi*200;OmegaP2=

5、2*pi*400;%带通截止频率OmegaS1=2*pi*100;OmegaS2=2*pi*600;%1.5kHz阻带Rp=2;%波纹系数As=20;%阻带衰减dBFs=2*103;%抽样频率2khzOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs;ws=OmegaS/Fs;%等效数字频率OmegaP_t=2*Fs*tan(wp/2);OmegaS_t=2*Fs*tan(ws/2);N,OmegaC=cheb1ord(OmegaP_t,OmegaS_t,Rp,As,s) b,a=cheby1(N,Rp,OmegaC,s);%A

6、F系统函数的分子 分母bz,az=bilinear(b,a,Fs)%双线性变换法 AF to DFw0=wp,ws;%四个频点Hx=freqz(bz,az,w0);%计算四个频点上对应的幅度响应H,w=freqz(bz,az);%计算0pi上的响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化并求dBdbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应plot(w/pi,dbH);%画图xlabel(数字频率域频率Omega/pi);ylabel(dB);axis(0,1,-300,5);grid运行,得N = 2;OmegaC =1.0

7、e+03 *1.2997 2.9062w0 = 0.6283 1.2566 0.3142 1.8850bz = 0.0512 -0.0000 -0.1024 -0.0000 0.0512az = 1.0000 -2.0733 2.4881 -1.5944 0.6125dbHx = 1.9997 1.9997 24.1325 22.3051图 2 双线性变换法设计IIR带通滤波器由程序返回得到的数值可以得知，这是一个2阶带通滤波器。转换为模拟原型滤波器的频率见OmegaC返回的值。数字滤波器的边界频率见w0返回的值。由bz,az返回值可得由dbHx返回值可知在四个边界频率处幅度响应大小:对应到模

8、拟频率，即有2004001006001.99971.999724.132522.3051满足题目要求的设计指标，通带波纹不大于2dB, 截止频率处衰减不小于20dB。7.14设计一个数字切比雪夫I型带阻滤波器，给定指标为：(1) 衰减,当 (2) 波纹，当 ， (3) 抽样频率试用双线性变换法进行设计，最后写出的表达式，并画出系统的幅频响应特性。解：双线性变换法设计程序如下%冲激响应不变法，ex714.m%设计数字切比雪夫I型带阻滤波器%双线性变换法clc;clear allOmegaP1=2*pi*500;OmegaP2=2*pi*3000;%带通截止频率OmegaS1=2*pi*1000;

9、OmegaS2=2*pi*2000;%阻带截止频率Rp=3;%通带波纹dBAs=30;%阻带衰减dBFs=10*103;%抽样频率10khzOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs;ws=OmegaS/Fs;%等效数字频率OmegaP_t=2*Fs*tan(wp/2);OmegaS_t=2*Fs*tan(ws/2);N,OmegaC=cheb1ord(OmegaP_t,OmegaS_t,Rp,As,s)%AF阶数和截至频率b,a=cheby1(N,Rp,OmegaC,stop,s);%AF系统函数的分子 分母bz,az

10、=bilinear(b,a,Fs)%双线性变换法F to DFw0=wp,ws;%四个频点Hx=freqz(bz,az,w0);%计算两个频点上对应的幅度响应H,w=freqz(bz,az);%计算0pi上的响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化并求dBdbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应plot(w/pi,dbH);%画图xlabel(数字频率域频率Omega/pi);ylabel(dB);axis(0,1,-150,5);grid运行,得N = 3;OmegaC = 1.0e+04 * 0.3430 2.7

11、528w0 = 0.3142 1.8850 0.6283 1.2566bz = 0.0946 -0.3508 0.7174 -0.8802 0.7174 -0.3508 0.0946az = 1.0000 -1.4600 0.3178 -0.3506 0.6885 0.2289 -0.3824dbHx = 0.2548 3.0000 39.8925 39.8925图 3双线性变换法设计IIR带阻滤波器由程序返回得到的数值可以得知，这是一个2阶带阻滤波器。转换为模拟原型滤波器的频率见OmegaC返回的值。数字滤波器的边界频率见w0返回的值。由bz,az返回值可得由dbHx返回值可知在四个边界频率

12、处幅度响应大小:对应到模拟频率，即有5003000100020000.25483.000039.892539.8925满足题目要求的设计指标，通带波纹不大于3dB, 截止频率处衰减大于30dB。7.15设计一个数字切比雪夫型带阻滤波器，给定指标为：(4) 衰减,当 (5) 波纹，当 ， (6) 抽样频率试用双线性变换法进行设计，最后写出的表达式，并画出系统的幅频响应特性。设计时请先想一想，这一题和上一题有什么相似处。由此应该得出什么结论。解：分析：这一题的各频率指标要求均是7.14的10倍，抽样频率也是10倍。此时对应的数字频率是一样的。如果用切比雪夫型滤波器，双线性变换法设计出来的与7.14

13、中的肯定是一样的（见图4）.这说明只要给定的各数字频率参数一样，用同种方法设计出来的数字滤波器是一样的。程序如下：%设计数字切比雪夫型带阻滤波器%双线性变换法clc;clear allOmegaP1=2*pi*5000;OmegaP2=2*pi*30000;%带通截止频率OmegaS1=2*pi*10000;OmegaS2=2*pi*20000;%阻带截止频率Rp=3;%通带波纹dBAs=30;%阻带衰减dBFs=100*103;%抽样频率10khzOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs;ws=OmegaS/Fs;%

14、等效数字频率OmegaP_t=2*Fs*tan(wp/2);OmegaS_t=2*Fs*tan(ws/2);N,OmegaC=cheb2ord(OmegaP_t,OmegaS_t,Rp,As,s)%AF阶数和截至频率b,a=cheby2(N,As,OmegaS,stop,s);%AF系统函数的分子 分母bz,az=bilinear(b,a,Fs)%双线性变换法F to DFw0=wp,ws;%四个频点Hx=freqz(bz,az,w0);%计算两个频点上对应的幅度响应H,w=freqz(bz,az);%计算0pi上的响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化

15、并求dBdbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应plot(w/pi,dbH);%画图xlabel(数字频率域频率Omega/pi);ylabel(dB);axis(0,1,-150,5);grid图 4双线性变换法设计IIR带阻滤波器（fs=100kHz）7.17要求设计一个数字带通滤波器，其抽样频率，通达截止频率为，通带衰减，阻带截止频率为，阻带衰减为。(1) 利用MATLAB工具箱中ellipord及ellip设计椭圆函数滤波器；(2) 利用MATLAB工具箱中cheblord及cheb1设计切比雪夫型滤波器；(3) 利用MATLAB工具箱中cheb

16、2ord及cheb2设计切比雪夫型滤波器；(4) 利用buttord及butter设计巴特沃思型滤波器。要求每种设计都给出系统函数并画出幅频特性（dB）、相频特性以及单位冲激响应。解：(1) 利用MATLAB工具箱中ellipord及ellip设计椭圆函数滤波器：核心语句为：%设计椭圆带通滤波器%双线性变换法clc;clear allOmegaP1=2*pi*5000;OmegaP2=2*pi*7000;%带通截止频率OmegaS1=2*pi*3500;OmegaS2=2*pi*8500;%阻带截止频率Rp=0.5;%通带波纹dBAs=45;%阻带衰减dBFs=35*103;%抽样频率35kh

17、zOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs/pi;ws=OmegaS/Fs/pi;%等效数字频率w0=wp,ws;%四个频点N,wc=ellipord(wp,ws,Rp,As)%DF阶数和截至频率bz,az=ellip(N,Rp,As,wc)%DF系统函数的分子 分母Hx=freqz(bz,az,w0*pi);%计算四个频点上对应的幅度响应H,w=freqz(bz,az);%计算01上的幅频响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化并求dBdbH=20*log10(abs(H)/ma

18、x(abs(H);%归一化的频率响应pha=unwrap(angle(H);%计算01上的相频响应figure(1)plot(w/pi,dbH);%画幅频图xlabel(数字频率域频率Omega/pi);ylabel(dB);axis(0,1,-80,5);gridfigure(2)plot(w/pi,pha);%画相频图xlabel(数字频率域频率Omega/pi);ylabel(Phasedegree);% axis(0,1,-150,5);griddisp(系统传递函数H(z);printsys(bz,az,z);figure(3)h=dimpulse(bz,az);stem(h)xla

19、bel(n);ylabel(Impulse response);图 5直接设计法设计IIR椭圆带通滤波器幅度响应图 6直接设计法设计IIR椭圆带通滤波器相位响应图 7直接设计法设计IIR椭圆带通滤波器单位冲激响应得到四个截止频率处的幅度响应如下：3.5578.50.20000.28570.40000.48570.49990.400047.738245.0760满足题目设定的指标要求。(2) 利用MATLAB工具箱中cheb1ord及cheby1设计切比雪夫型滤波器:%设计chebshev I带通滤波器clc;clear allOmegaP1=2*pi*5000;OmegaP2=2*pi*700

20、0;%带通截止频率OmegaS1=2*pi*3500;OmegaS2=2*pi*8500;%阻带截止频率Rp=0.5;%通带波纹dBAs=45;%阻带衰减dBFs=35*103;%抽样频率35khzOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs/pi;ws=OmegaS/Fs/pi;%等效数字频率w0=wp,ws;%四个频点N,wc=cheb1ord(wp,ws,Rp,As)%DF阶数和截至频率bz,az=cheby1(N,Rp,wc);%DF系统函数的分子 分母Hx=freqz(bz,az,w0*pi);%计算四个频点上对

21、应的幅度响应H,w=freqz(bz,az);%计算01上的幅频响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化并求dBdbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应pha=unwrap(angle(H);%计算01上的相频响应figure(1)plot(w/pi,dbH);%画幅频图xlabel(数字频率域频率Omega/pi);ylabel(dB);axis(0,1,-150,5);gridfigure(2)plot(w/pi,pha);%画相频图xlabel(数字频率域频率Omega/pi);ylabel(Phasedegr

22、ee);% axis(0,1,-150,5);griddisp(系统传递函数H(z);printsys(bz,az,z);figure(3)h=dimpulse(bz,az);stem(h)xlabel(n);ylabel(Impulse response);grid 图 8直接设计法设计IIR切比雪夫型带通滤波器幅度响应图 9直接设计法设计IIR切比雪夫型带通滤波器相位响应图 10直接设计法设计IIR切比雪夫型带通滤波单位冲激响应得到四个截止频率处的幅度响应如下：3.5578.50.20000.28570.40000.48570.49990.499962.451351.1270满足题目设定的

23、指标要求。(3) 利用MATLAB工具箱中cheb2ord及cheby2设计切比雪夫型滤波器：核心语句为：%设计chebshev 带通滤波器clc;clear allOmegaP1=2*pi*5000;OmegaP2=2*pi*7000;%带通截止频率OmegaS1=2*pi*3500;OmegaS2=2*pi*8500;%阻带截止频率Rp=0.5;%通带波纹dBAs=45;%阻带衰减dBFs=35*103;%抽样频率35khzOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs/pi;ws=OmegaS/Fs/pi;%等效数字频

24、率N,wc=cheb2ord(wp,ws,Rp,As)%AF阶数和截至频率bz,az=cheby2(N,As,wc);%AF系统函数的分子 分母w0=wp,ws;%四个频点Hx=freqz(bz,az,w0*pi);%计算两个频点上对应的幅度响应H,w=freqz(bz,az);%计算01上的幅频响应dbHx=-20*log10(abs(Hx)/max(abs(H)%归一化并求dBdbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应pha=unwrap(angle(H);%计算01上的相频响应figure(1)plot(w/pi,dbH);%画幅频图xlabel(数

25、字频率域频率Omega/pi);ylabel(dB);axis(0,1,-100,5);gridfigure(2)plot(w/pi,pha);%画相频图xlabel(数字频率域频率Omega/pi);ylabel(Phasedegree);% axis(0,1,-150,5);griddisp(系统传递函数H(z);printsys(bz,az,z);figure(3)h=dimpulse(bz,az);stem(h)xlabel(n);ylabel(Impulse response);grid图 11直接设计法设计IIR切比雪夫型带通滤波器幅度响应图 12 直接设计法设计IIR切比雪夫型带

26、通滤波器相位响应图 Error! Main Document Only.直接设计法设计IIR切比雪夫型带通滤波器单位冲激响应得到四个截止频率处的幅度响应如下：3.5578.50.20000.28570.40000.48570.00090.530245.000045.0000满足题目设定的指标要求。(4) 利用buttord及butter设计巴特沃思型滤波器。核心语句为：%设计chebshev 带通滤波器clc;clear allOmegaP1=2*pi*5000;OmegaP2=2*pi*7000;%带通截止频率OmegaS1=2*pi*3500;OmegaS2=2*pi*8500;%阻带截止

27、频率Rp=0.5;%通带波纹dBAs=45;%阻带衰减dBFs=35*103;%抽样频率35khzOmegaP=OmegaP1,OmegaP2;OmegaS=OmegaS1,OmegaS2;wp=OmegaP/Fs/pi;ws=OmegaS/Fs/pi;%等效数字频率w0=wp,ws;%四个频点N,wc=buttord(wp,ws,Rp,As)%AF阶数和截至频率bz,az=butter(N,wc);%AF系统函数的分子 分母Hx=freqz(bz,az,w0*pi);%计算两个频点上对应的幅度响应H,w=freqz(bz,az);%计算01上的幅频响应dbHx=-20*log10(abs(H

28、x)/max(abs(H)%归一化并求dBdbH=20*log10(abs(H)/max(abs(H);%归一化的频率响应pha=unwrap(angle(H);%计算01上的相频响应figure(1)plot(w/pi,dbH);%画幅频图xlabel(数字频率域频率Omega/pi);ylabel(dB);axis(0,1,-200,5);gridfigure(2)plot(w/pi,pha);%画相频图xlabel(数字频率域频率Omega/pi);ylabel(Phasedegree);% axis(0,1,-150,5);griddisp(系统传递函数H(z);printsys(bz

29、,az,z);figure(3)h=dimpulse(bz,az);stem(h)xlabel(n);ylabel(Impulse response);grid图 Error! Main Document Only.直接设计法设计IIR巴特沃思带通滤波器幅度响应图 15直接设计法设计IIR巴特沃思带通滤波器相位响应图 16直接设计法设计IIR巴特沃思带通滤波器单位冲激响应得到四个截止频率处的幅度响应如下：3.5578.50.20000.28570.40000.48570.10540.105461.833145.0000满足题目设定的指标要求。FIR数字滤波器设计作业周昂 通信工程 201304

30、0150278.2设计一个线性相位FIR数字高通滤波器。技术指标为，。画出及及曲线。解：分析：由于要求阻带最小衰减为55dB, 故考虑用blackman窗进行设计。程序如下：%8.2程序，FIR高通滤波器clc,clear all;wp=0.7*pi;ws=0.5*pi;dw=abs(ws-wp);%过渡带宽N0=ceil(11*pi/dw);%计算blackman窗长N=N0+mod(N0+1,2);%确保N是奇数n=0:N-1;wd=(blackman(N);%求blakman函数wc=(wp+ws)/2/pi;%理想低通滤波器的截止频率hn=fir1(N-1,wc,high,blackm

31、an(N);H,w=freqz(hn,1,1000);%求频率响应dbH=20*log10(abs(H);dw=2*pi/1000;As=-max(dbH(1:wp/dw+1)%检查通带最大衰减Rp=-min(dbH(ws/dw+1:501)%检查阻带最小衰减pha=unwrap(angle(H);%计算01上的相频响应figure(1);stem(n,hn);title(单位冲激响应(db);xlabel(n);ylabel(h(n);gridfigure(2)subplot(2,2,1 3)plot(w/pi,dbH);title(幅度响应(db);xlabel(omage/pi);yla

32、bel(20log|H(ej/omega)|(dB);axis(0,1,-200,5);grid;subplot(2,2,2)plot(w/pi,dbH);title(幅度响应(db);xlabel(omage/pi);ylabel(20log|H(ej/omega)|(dB);axis(0.5,0.9,-10,10);grid;set(gca,xtickmode,manual,xtick,0.5:0.1:0.9);set(gca,ytickmode,manual,ytick,-10,-3,0,3,10);subplot(2,2,4)plot(w/pi,dbH);title(幅度响应(db);

33、xlabel(omage/pi);ylabel(20log|H(ej/omega)|(dB);grid;axis(0.3 0.6 -100 -50)set(gca,xtickmode,manual,xtick,0.3:0.1:0.6);set(gca,ytickmode,manual,ytick,-100,-80,-55,-50);figure(3);plot(w/pi,pha);%画幅频图xlabel(数字频率域频率Omega/pi);ylabel(Phase/degree);% axis(0,1,-150,5);grid运行,得As =72.3744Rp = 0.0023图 17 加bla

34、ckman窗的FIR高通滤波器单位冲激响应图 18 加blakcman窗的FIR高通滤波器幅频响应图 19加blackman窗的FIR高通滤波器相频响应由运行结果可以看出，阻带最小衰减达到了72.37，满足设计指标。由幅频图可以看出，高通滤波器在取最小衰减，已经达到了题目设计的55dB的指标。由相频响应可以看出，在通带满足线性相位。信号得到无失真传输。8.3设计一个线性相位FIR数字带通滤波器。技术指标为，。画出及及曲线。解：分析：由于要求阻带最小衰减为75dB,故考虑用kaiser窗进行设计。用kaiser窗设计时需要根据阻带最小衰减值As计算出beta值，这时kaiser窗与其他窗在使用上

35、的一点区别。程序如下：%8.3程序，FIR带通滤波器clc,clear all;wp1=0.4*pi;wp2=0.5*pi;ws1=0.2*pi;ws2=0.7*pi;As0=75;B=min(wp1-ws1,ws2-wp2);%过渡带宽N0=ceil(As0-7.95)/(2.286*B)+1;%计算kaiser窗长N=N0+mod(N0+1,2);%确保N是奇数n=0:N-1;beta=0.1102*(As0-8.7);%betawd=(kaiser(N,beta);%求kaiser函数wc1=(wp1+ws1)/2/pi;wc2=(wp2+ws2)/2/pi;wc=wc1,wc2;hn=

36、fir1(N-1,wc,bandpass,kaiser(N,beta);H,w=freqz(hn,1,1000);%求频率响应dbH=20*log10(abs(H);dw=pi/1000;Rp=-min(dbH(wp1/dw+1:1:wp2/dw+1)%检查通带最大衰减As=-max(max(dbH(1:1:ws1/dw+1),max(dbH(ws2/dw+1:1:1000)%检查阻带最小衰减pha=unwrap(angle(H);%计算01上的相频响应figure(1);stem(n,hn);%画冲激响应图title(单位冲激响应(db);xlabel(n);ylabel(h(n);grid

37、figure(2)%画幅频图subplot(2,2,1 3)plot(w/pi,dbH);xlabel(Omega/pi);ylabel(20log|H(ej/omega)|(dB);grid;subplot(2,2,2)plot(w/pi,dbH);xlabel(Omega/pi);ylabel(20log|H(ej/omega)|(dB);grid;axis(0.3 0.6 -10 10)set(gca,xtickmode,manual,xtick,0.3,0.4,0.5,0.6);set(gca,ytickmode,manual,ytick,-10,-3,0,3,10);subplot(

38、2,2,4)plot(w/pi,dbH);xlabel(Omega/pi);ylabel(20log|H(ej/omega)|(dB);grid;axis(0.6 0.9 -100 -50)set(gca,xtickmode,manual,xtick,0.6,0.7,0.8,0.9);set(gca,ytickmode,manual,ytick,-100,-75,-50);figure(3);plot(w/pi,pha);%画相频图xlabel(Omega/pi);ylabel(Phase/degree);grid运行,得As =72.8692Rp = 0.0021图 20 加blackman

39、窗的FIR带通滤波器单位冲激响应图 21 加blackman窗的FIR带通滤波器幅频响应图 22加blackman窗的FIR带通滤波器相频响应可以看出，设定阻带衰减为75dB时，并不满足要求。因此将衰减增至80dB,得到返回的阻带最大衰减为75.9251dB，满足要求。由相频响应可以看出，在通带满足线性相位。信号得到无失真传输。8.4设计一个线性相位FIR数字带阻滤波器。技术指标为，。画出及及曲线。解：分析：由于要求阻带最小衰减为80dB,故考虑用kaiser窗进行设计。程序如下：%8.4程序，FIR带阻滤波器%通带截止频率wp 0.35pi 0.8pi%阻带截止频率ws 0.5pi 0.65pi 最小衰减As80dBclc,clear all;wp1=0.35*pi;wp2=0.8*pi;ws1=0.5*pi;ws2=0.65*pi;As0=80;B=min(ws1-wp1,wp2-ws2);%过渡带宽N0=ceil(As0-7.95)/