数字滤波器的设计步骤及程序实现湖南理工学院信息与通信工程学院一、IIR 脉冲响应不变法设计步骤1、已知实际数字指标as s ap p ，，，ωω2、将数字指标化为原型模拟指标As s Ap p ，，，ΩΩ，可设T=pi, T /ω=Ω3、求原型模拟滤波器的c N Ω，，其中：⎥⎥⎤⎢⎢⎡ΩΩ--=)/lg(2)]110/()110lg[(10/10/sp A A s p NNA pcp p 210/110-Ω=Ω NA scs s 210/110-Ω=Ω ][cs cp c ΩΩ∈Ω，4、根据N 写出归一化原型系统函数)(p G a5、用c s p Ω=/代入得原型系统函数cs p a a p G s H Ω==/)()(6、将)(s H a 化为部分分式展开形式∑-=kka s s A s H )(7、写出)(z H 的极点Ts k k ez =，并写出)(z H 的部分分式展开形式∑--⋅=11)(z z A T z H kk8、将)(z H 化为分子分母形式，验证设计结果。

二、IIR 双线性变换法设计步骤1、已知实际数字指标as s ap p ，，，ωω2、将数字指标化为原型模拟指标As s Ap p ，，，ΩΩ，可设T=2, 2tan 2ω⋅=ΩT 3、求原型模拟滤波器的c N Ω，，其中：⎥⎥⎤⎢⎢⎡ΩΩ--=)/lg(2)]110/()110lg[(10/10/s p A A s p NNA pcp p 210/110-Ω=Ω NA scs s 210/110-Ω=Ω ][cs cp c ΩΩ∈Ω，4、根据N 写出归一化原型系统函数)(p G a5、用c s p Ω=/代入得原型系统函数cs p a a p G s H Ω==/)()(6、用11112--+-⋅=Z Z T s 代入原型系统函数)(s H a 得11112)()(--+-⋅==Z Z Ts a s H z H 8、将)(z H 整理成分子分母形式，验证设计结果。

三、FIR 窗函数法设计步骤1、已知实际数字指标as s ap p ，，，ωω2、根据as 选窗的类型：矩形窗as<21dB, A=1.8π,窗函数是 boxcar(N)；三角窗as<25dB, A=6.1π,窗函数是 bartlett(N)；汉宁窗as<44dB, A=6.2π,窗函数是 hanning(N)；哈明窗as<53dB, A=6.6π,窗函数是 hamming(N)；布莱克曼窗as<74dB, A=11π,窗函数是 blackman(N)。

3、根据过渡带p s B t ωω-=和窗类型求总点数t B A N /≈。

4、根据2/)(s p c ωωω+=写出理想频响指标)()()(ωθωωj dg j d e H e H ⋅=5、根据)(ωj d eH 算出ωπωππωd e e H n h n j j d d ⎰-⋅=)(21)( 6、对)(n h d 加窗得设计结果)()()(n w n h n h d ⋅=8、写出∑-=nzn h z H )()(，验证设计结果。

四、FIR 频率采样法设计步骤1、已知实际数字指标as s ap p ，，，ωω2、根据as 选过渡带点数m3、根据过渡带p s B t ωω-=和过渡带点数m 求总点数t B m N /2)1(π⋅+≥。

4、根据c ω求出πω2Nk c c ⋅=，设置过渡值)1 0(，∈T 5、根据约束条件构建理想频响的采样指标)()()(k j g d e k H k H θ⋅=6、对)(k H d 进行IDFT 变换得)(n h ，取实部。

7、写出∑-=nzn h z H )()(，验证设计结果，优化过渡值大小、过渡点位置和过渡点多少。

一、IIR滤波器设计：脉冲响应不变法实现程序%用脉冲响应不变法设计butterworth数字低通滤波器%技术指标：wp=0.3*pi rad, ap=2dB, ws=0.5*pi rad, as=10dBclc; clear; close all; format compact;%程序初始化wp=0.3*pi, ap=2, ws=0.5*pi, as=10,%输入数字指标T=pi,%假设采样周期，用于设计原型模拟滤波器，不影响H(z)的设计结果Wp=wp/T, Ap=ap, Ws=ws/T, As=as,%将数字指标转化为原型模拟指标M=log10( (10 .^ (0.1*Ap) - 1)./(10 .^ (0.1*As) - 1) ) / ...(2*log10(Wp/Ws)) ,%计算滤波器阶数N = ceil( M),%滤波器阶数向上取整Wcp = Wp / ( (10^(.1*Ap) - 1)^(1/(2*N))),%通带边界精确满足的截止频率Wcs = Ws / ( (10^(.1*As) - 1)^(1/(2*N))),%阻带边界精确满足的截止频率Wc=Wcp,%截止频率用通带边界精确满足的截止频率%Wc=(Wcp+Wcs)/2,%通带阻带边界都有余量的截止频率%Wc=Wcs,%截止频率用阻带边界精确满足的截止频率[bp,ap]=butter(N,1,'s'),%求归一化原型滤波器系统函数Ga(p)P157tf(bp,ap,'variable','p'),%显示Ga(p)[bs,as]=lp2lp(bp,ap,Wc),%去归一化得原型滤波器系统函数Ha(s)tf(bs,as),%显示Ha(s)，分子不足前面补0[Ak,sk]=residue(bs,as),%将Ha(s)按部分分式形式展开ak=T*Ak,zk=exp(sk*T),%将Ha(s)的部分分式参数转换为H(z)的部分分式参数[bz,az]=residuez(ak,zk,0),%将H(z)的部分分式形式化为分子分母等阶形式tf(bz,az,'variable','z^-1'),%显示系统函数%[bz1,az1] = impinvar(bs,as,1/T)%调用impinvar函数验证%tf(bz1,az1,'variable','z^-1'),%显示验证系统函数,z^-1式的分子不足是后面补0freqz(bz,az,100),%绘出频率特性曲线，检验设计指标二、IIR滤波器设计：双线性变换法实现程序%用双线性变换法设计butterworth数字低通滤波器%技术指标：wp=0.3*pi rad, ap=2dB, ws=0.5*pi rad, as=10dBclc; clear; close all; format compact;%程序初始化wp=0.3*pi, ap=2, ws=0.5*pi, as=10,%输入数字指标T=2,%假设采样周期，用于设计原型模拟滤波器，不影响H(z)的设计结果Wp=(2/T)*tan(wp/2), Ap=ap,Ws=(2/T)*tan(ws/2),As=as,%将数字指标预畸变成原型模拟指标M=log10( (10 .^ (0.1*Ap) - 1)./(10 .^ (0.1*As) - 1) ) / ...(2*log10(Wp/Ws)) ,%计算滤波器阶数N = ceil( M),%滤波器阶数向上取整Wcp = Wp / ( (10^(.1*Ap) - 1)^(1/(2*N))),%通带边界精确满足的截止频率Wcs = Ws / ( (10^(.1*As) - 1)^(1/(2*N))),%阻带边界精确满足的截止频率Wc=Wcp,%截止频率用通带边界精确满足的截止频率%Wc=(Wcp+Wcs)/2,%通带阻带边界都有余量的截止频率%Wc=Wcs,%截止频率用阻带边界精确满足的截止频率[bp,ap]=butter(N,1,'s'),%求归一化原型滤波器系统函数Ga(p)P157tf(bp,ap,'variable','p'),%显示Ga(p)[bs,as]=lp2lp(bp,ap,Wc),%去归一化得原型滤波器系统函数Ha(s)tf(bs,as),%显示Ha(s)，分子不足前面补0[bz,az] = bilinear (bs,as,1/T),%将模拟低通原型转换为数字低通tf(bz,az,'variable','z^-1'),%显示系统函数freqz(bz,az,100),%绘出频率特性曲线，检验设计指标三、FIR滤波器设计：窗函数法实现程序%用窗函数法设计数字低通滤波器%技术指标：wp=0.27*pi rad, ap=2dB, ws=0.40*pi rad, as=10dB。

clc; clear; close all; format compact;%程序初始化wp=0.27*pi, ap=2, ws=0.40*pi, as=10,%输入数字指标%根据as选择窗函数的类型并输入参数A，计算窗口长度M%矩形窗as<21dB,A=1.8*pi,窗函数是boxcar(N)%三角窗as<25dB,A=6.1*pi,窗函数是bartlett(N)%汉宁窗as<44dB,A=6.2*pi,窗函数是hanning(N)%哈明窗as<53dB,A=6.6*pi,窗函数是hamming(N)%布莱克曼窗as<74dB,A=11*pi,窗函数是blackman(N)A=1.8*pi,%因as=10dB选矩形窗Bt=ws-wp; %计算过渡带宽M=ceil(A/Bt);%根据窗函数的类型计算长度if mod(M,2)==0; N=M+1, else N=M, end; %选用第一类滤波器wc=(wp+ws)/2, %转折频率一般取通带频率和阻带频率的中点n=-30:40;r=(N-1)/2; %用于计算理想低通单位脉冲响应中数据hdn=sin(wc*((n-r)+eps))./(pi*((n-r)+eps)); %参见教材P202求理想单位脉冲响应wn=boxcar(N); %窗函数数据m=0:N-1;r=(N-1)/2; hm=sin(wc*((m-r)+eps))./(pi*((m-r)+eps));hn=hm'.*wn; %理想单位脉冲响应加窗处理figure(1),%绘加窗处理过程图subplot(3,1,1),stem(n,hdn,'r.'),grid on ,axis([-15,30,-0.2,0.5])subplot(3,1,2),stem([0:N-1],wn,'.'),grid on ,axis([-15,30,-0.4,1.4])subplot(3,1,3),stem([0:N-1],hn,'k.'),grid on ,axis([-15,30,-0.2,0.5])figure(2);freqz(hn,1,100);% 绘频率特性曲线图，检验设计指标figure(3),%绘幅度响应函数Hg(ω)图Hejw=fft(hn,256); %计算频率响应函数k=[0:255];wk=2*pi/256*k;Hgw=(Hejw.').*exp(j*wk*(N-1)/2); %计算幅频响应plot(wk/pi,real(Hgw));xlabel('ω/π');ylabel('Hg(ω)');%绘图四、FIR滤波器设计：频率采样法实现程序%用频率采样法设计FIR低通滤波器%技术指标：wc=0.3*pi rad, N=15，不加过渡点clc; clear; close all; format compact;%程序初始化wc=0.3*pi;N=15;%设计阶数N为奇的第一类滤波器%根据约束条件确定H(k)的值k=[0:N-1],w=2*pi/N*k;kc=fix(wc*N/(2*pi)),%求频点及转折频率对应的k值Hk_abs=[ones(1,kc+1),zeros(1,N-2*kc-1),ones(1,kc)],%采样频点幅值Hk_angles=-(N-1)/N*pi*k;%采样频点相位Hk= Hk_abs.*exp(j*Hk_angles);%采样频点的H(k)hn=real(ifft(Hk)),%求H(k)的IDFT得单位脉冲响—即设计结果%理解频域离散与时域的周期延拓n=-30:40;r=(N-1)/2;hdn1=sin(wc*((n-r)+eps))./(pi*((n-r)+eps)); %参见教材P202hdn2=sin(wc*((n-r)+N+eps))./(pi*((n-r)+N+eps)); %参见教材P202hdn3=sin(wc*((n-r)-N+eps))./(pi*((n-r)-N+eps)); %参见教材P202figure(1),%绘时域的周期延拓叠加图subplot(4,1,1),stem(n,hdn1,'r.'),grid on ,axis([-15,30,-0.1,0.4])subplot(4,1,2),stem(n,hdn2,'.'),grid on ,axis([-15,30,-0.1,0.4])subplot(4,1,3),stem(n,hdn3,'.'),grid on ,axis([-15,30,-0.1,0.4])subplot(4,1,4),stem([0:14],hn,'k.'),grid on ,axis([-15,30,-0.1,0.4])figure(2)%绘幅度谱变化过程图w1=[0,wc,wc+eps,2*pi-wc,2*pi-wc+eps,2*pi]/pi;xk=[1,1,0,0,1,1]plot(w1,xk,'r:'),xlabel('ω/π'),axis([0,2,-0.3,1.2])%绘理想幅频曲线hold on,stem(2/N*k, Hk_abs,'.'),%绘频率采样图%绘设计结果的幅频响应图Hejw=fft(hn,256); %计算频率响应K=[0:255];wk=2*pi/256*K;Hgw=Hejw.*exp(j*wk*(N-1)/2); %计算幅频响应hold on;plot(wk/pi,real(Hgw),'k');xlabel('ω/π');ylabel('Hg(ω)');%绘图。