《数字信号处理》课程研究性学习报告 试点班专用 姓名 学号 同组成员 指导教师 时间 数字滤波器设计专题研讨 【目的】 (1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。

(2) 掌握各种窗函数的时频特性及对滤波器设计的影响。

(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。

【研讨题目】 基本题 1．IIR 数字滤波器设计 设计一个IIR数字低通滤波器，其能取代下列指标的模拟低通滤波器(系统的抽样频率为44.1kHz) fp=2kHz， fs=10kHz , Ap=0.5dB, As=50dB (1) 分别用双线性变换和冲激响应不变法设计一个BW型数字低通滤波器，并进行比较。

(2) 用双线性变换分别设计Chebyshev I型Chebyshev I I型和椭圆型数字低通滤波器，并进行比较。

【温磬提示】 在数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。

模拟滤波器设计指标 【设计步骤】 数字低通滤波器 模拟低通滤波器 数字滤波器设计指标 频率变换 → 双线性法 冲击不变法 【仿真结果】 所设计滤波器的幅度响应和相位响应 BW型、Chebyshev I型、Chebyshev I I型和椭圆型滤波器的零极点分布 【结果分析】 双线性变换和冲激响应不变法所设计的滤波器的性能有什么不同。

BW型、Chebyshev I型、Chebyshev I I型和椭圆型滤波器的零极点分布各有什么特点。

答：
双线性法得到的模拟频率与数字频率是非线性的，但是消除了频谱混叠的误差 脉冲响应不变法的模拟频率域数字频率是线性的，但是有频谱混叠误差。

在极点图中，BW型极点离单位圆最远，椭圆极点离单位圆最近。因而BW的稳定性最好，椭圆的稳定性最差。

【自主学习内容】 极点图的做法 【阅读文献】 《信号与系统》 【发现问题】 (专题研讨或相关知识点学习中发现的问题)：
Cheby2型做的幅度响应在ws之后没有波动。

【仿真程序】 wp=2*pi*2000;ws=2*pi*10000;Ap=0.5;As=50;Fs=44100; Wp=wp/Fs;Ws=ws/Fs; wp1=2*Fs*tan(Wp/2);ws1=2*Fs*tan(Ws/2); [N,wc]=buttord(wp1,ws1,Ap,As,＇s＇); [num,den]=butter(N,wc,＇s＇); [numd,dend]=bilinear(num,den,Fs); w=linspace(0,pi,512); h=freqz(numd,dend,w); norm=max(abs(h)); numd=numd/norm; A=20*log10(abs(h)/norm); subplot(1,3,1);plot(w/pi,abs(h));xlabel(＇normalized frequency＇);ylabel(＇幅度响应＇);title(＇BW型 双线性法＇); subplot(1,3,2);plot(w/pi,unwrap(angle(h)));xlabel(＇normalized frequency＇);ylabel(＇相位响应＇);title(＇BW型 双线性法＇); subplot(1,3,3);plot(w/pi,A);xlabel(＇normalized frequency＇);ylabel(＇Gain in dB＇); title(＇BW型 双线性法＇); w=[Wp Ws]; h=freqz(numd,dend,w); A=20*log10(abs(h)/norm) 冲击响应不变法：
wp=2*pi*2000;ws=2*pi*10000;Ap=0.5;As=50;Fs=44100; Wp=wp/Fs;Ws=ws/Fs; [N,wc]=buttord(wp,ws,Ap,As,＇s＇); [num,den]=butter(N,wc,＇s＇); [numd,dend]=impinvar(num,den,Fs); w=linspace(0,pi,512); h=freqz(numd,dend,w); norm=max(abs(h)); numd=numd/norm; A=20*log10(abs(h)/norm); subplot(1,3,1);plot(w/pi,abs(h));xlabel(＇normalized frequency＇);ylabel(＇幅度响应＇);title(＇BW型 脉冲不变法＇); subplot(1,3,2);plot(w/pi,unwrap(angle(h)));xlabel(＇normalized frequency＇);ylabel(＇相位响应＇);title(＇BW型 脉冲不变法＇); subplot(1,3,3);plot(w/pi,A);xlabel(＇normalized frequency＇);ylabel(＇Gain in dB＇); title(＇BW型 脉冲不变法＇); w=[Wp Ws]; h=freqz(numd,dend,w); A=20*log10(abs(h)/norm) 【研讨题目】 基本题 2．窗函数研究 分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性，并进行比较。

【题目分析】 这几种窗函数的时域定义如下：
1. 矩形窗：
2. Hanning窗：
3. Hamming窗：
4. Blankman窗：
5. Kaiser窗：
取定N值，再做L点DFT就能得到它们的频域特性。

【仿真结果】 【结果分析】 无论是对矩形窗、汉纳窗、哈明窗、布莱克曼窗的横向比较，还是对取不同β值的凯泽窗的纵向比较，都发现，窗函数通过牺牲主瓣宽度来降低频率泄露现象。

【问题探究】 在谱分析中如何选择窗函数，在滤波器设计中如何选择窗函数？ 要考虑旁瓣对频谱的影响，并且尽量用能量集中的主瓣。

【仿真程序】 1.%Rectangle Window,Hanning Window,Hamming Window,Blackman Window N=64;L=1024; k=0:1:N-1; wRec=ones(1,N); wHan=0.5-0.5*cos(2*pi*k/N); wHam=0.54-0.46*cos(2*pi*k/N); wBla=0.42-0.5*cos(2*pi*k/N)+0.08*cos(4*pi*k/N); WRec=fftshift(fft(wRec,L)); WHan=fftshift(fft(wHan,L)); WHam=fftshift(fft(wHam,L)); WBla=fftshift(fft(wBla,L)); w=(2/L)*(0:L-1)-1; subplot(2,2,1);plot(w,abs(WRec)/max(abs(WRec))); title(＇The Amplitude Response of Rectangle Window＇); xlabel(＇Frequency(πHz)＇);ylabel(＇Normalized Amplitude＇); subplot(2,2,2);plot(w,abs(WHan)/max(abs(WHan))); title(＇The Amplitude Response of Hanning Window＇); xlabel(＇Frequency(πHz)＇);ylabel(＇Normalized Amplitude＇); subplot(2,2,3);plot(w,abs(WHam)/max(abs(WHam))); title(＇The Amplitude Response of Hamming Window＇); xlabel(＇Frequency(πHz)＇);ylabel(＇Normalized Amplitude＇); subplot(2,2,4);plot(w,abs(WBla)/max(abs(WBla))); title(＇The Amplitude Response of Blankman Window＇); xlabel(＇ Frequency(πHz)＇);ylabel(＇Normalized Amplitude＇); 2.%Blackman Window旁瓣放大图 N=64;L=1024; k=0:1:N-1; wBla=0.42-0.5*cos(2*pi*k/N)+0.08*cos(4*pi*k/N); WBla=fftshift(fft(wBla,L)); w=(2/L)*(0:L-1)-1; plot(w,abs(WBla)/max(abs(WBla))); title(＇Blackman Window旁瓣放大图＇); set(gca,＇XTick＇,-0.2:0.01:0.2); set(gca,＇YTick＇,0:0.0001:0.0015); axis([-0.2,0.2,0,0.0015]); grid on; xlabel(＇Frequency(πHz)＇);ylabel(＇Normalized Amplitude＇); 3.%Kaiser Window N=64;L=1024; k=0:1:N-1; n=1; while n<=4; b=2*(n-1); wKai=(kaiser(N,b))＇; WKai=fftshift(fft(wKai,L)); w=(2/L)*(0:L-1)-1; figure(n);plot(w,abs(WKai)/max(abs(WKai))); title(＇Kaiser window(β=)＇); xlabel(＇Frequency(πHz)＇);ylabel(＇Normalized Amplitude＇); n=n+1; end 【研讨题目】 基本题 3． 窗函数法设计FIR 数字滤波器 (1)分别用Blackman窗和Kaiser窗法设计一个满足下列指标的线性相位的FIR低通滤波器 Wp=0.4p rad, Ap=0.5 dB, Ws=0.6p rad, As=55dB (2)（M5-5）在用窗口法设计FIR滤波器时，由于理想滤波器的幅度响应在截频处发生突变，使得设计出的滤波器的幅度响应发生振荡，这个现象被称为Gibbs现象。解决这个问题的一个方案是本书中介绍的用逐步衰减的窗函数。另一个方案是使理想滤波器过渡带为渐变的，如下图所示具有线性过渡带的理想低通滤波器的频率响应，试用窗口法设计逼近该频率响应的FIR滤波器。

题3图 【(2)单位脉冲响应证明】 试证该滤波器的单位脉冲响应为 其中：， 【设计步骤】 (1)根据给定的滤波器技术指标，选择滤波器的阶数N并分别应用两种窗函数w(n). 即用一个有限长度的窗函数序列w(n)来截取一个无限长的序列hd（n），获得一个 有限长序列h（n），即h（n）=w（n）*hd（n）。

(2) 渐变的窗选为hamming窗。取wp=0.55pi,ws=0.45pi,As=25db,Ap=1db。

设hamming窗的长度为M,矩形窗的长度为M1。M1取不同的值作图，其中M=7。

【仿真结果】 所设计滤波器的幅度响应和相位响应 【结果分析】 1. 由不同窗函数设计的滤波器的幅度响应我们可以看出，Blackman窗设计的滤波器在通阻带都有波动。且阻带衰减值固定。而Kaiser窗法能够根据不同的设计要求来改变窗函数的长度从而来调节通阻带的衰减，使最终的滤波器比较符合实际要求。

2. 通过逐步衰减的窗函数法和渐变过渡带两种方法设计的滤波器都能较好的设计较小Gibbs现象的滤波器。通过理论分析我们可以知道，FIR滤波器的波动是由于窗函数的主瓣面积和旁瓣面积决定的，所以增加采样点并不能改变通阻带的波动。而对于滤波器的渐变法，则可以通过改变点数很好的减少波动的范围。

【自主学习内容】 【阅读文献】 【发现问题】 (专题研讨或相关知识点学习中发现的问题)：
【问题探究】 通过实验讨论如何控制滤波器的阻带衰减 【仿真程序】 (1) 用Blackman窗 Ap=0.5;As=45;Wp=0.4*pi;Ws=0.6*pi; N=floor(11.4*pi/(Ws-Wp)); N=mod(N+1,2)+N;M=N-1;w=blackman(N)＇;Wc=(Wp+Ws)/2;k=0:M; hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);h=hd.*w; omega=linspace(0,pi,512);mag=freqz(h,[1],omega); figure plot(omega/pi,20*log10(abs(mag)));grid on; xlabel(＇Normalized frequency＇); ylabel(＇Gain in dB＇); title(＇Blackman＇) 用Kaiser窗 wp=0.4*pi;ws=0.6*pi;Ap=0.5;As=55;f=[wp/pi ws/pi];a=[1,0]; dev=10^(-As/20)*ones(1,length(a)); [M,Wc,beta,ftype]=kaiserord(f,a,dev);M=mod(M,2)+M;h=fir1(M,Wc,ftype,kaiser(M+1,beta)); omega=linspace(0,pi,512);mag=freqz(h,[1],omega); figure plot(omega/pi,20*log10(abs(mag)));grid on zoom on xlabel(＇Normalized frequency＇);ylabel(＇Gain in dB＇);title(＇ Kaiser＇) (2) wp=0.55*pi;ws=0.45*pi;Ap=1;As=25; N=ceil(7*pi/(wp-ws)); N=mod(N+1,2)+N;M=N-1;w=hamming(N);wc=(wp+ws)/2;k=0:M; hd=(wc/pi)*sinc(wc*(k-0.5*M)/pi);h=hd＇.*w; omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=abs(mag); plot(omega/pi,magdb,＇b＇);grid;w=ws-wp; M1=8;k2=-M1:M1;wc=(wp+ws)/2; hd=sinc(w*k2/2).*(sin(wc*k2)./(k2.*pi));hd(M1+1)=wc/pi; omega2=linspace(0,pi,512);mag2=freqz(hd,[1],omega2);magdb2=abs(mag2); hold on; plot(omega2/pi,magdb2,＇r＇); legend(＇逐步衰减＇,＇过渡带渐变＇); title(＇M1=8＇) wp=0.55*pi;ws=0.45*pi;Ap=1;As=25; N=ceil(7*pi/(wp-ws)); N=mod(N+1,2)+N;M=N-1;w=hamming(N);wc=(wp+ws)/2;k=0:M; hd=(wc/pi)*sinc(wc*(k-0.5*M)/pi);h=hd＇.*w; omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=abs(mag); plot(omega/pi,magdb,＇b＇);grid;w=ws-wp; M1=16;k2=-M1:M1;wc=(wp+ws)/2; hd=sinc(w*k2/2).*(sin(wc*k2)./(k2.*pi));hd(M1+1)=wc/pi; omega2=linspace(0,pi,512);mag2=freqz(hd,[1],omega2);magdb2=abs(mag2); hold on; plot(omega2/pi,magdb2,＇r＇); legend(＇逐步衰减＇,＇过渡带渐变＇); title(＇M1=16＇) wp=0.55*pi;ws=0.45*pi;Ap=1;As=25; N=ceil(7*pi/(wp-ws)); N=mod(N+1,2)+N;M=N-1;w=hamming(N);wc=(wp+ws)/2;k=0:M; hd=(wc/pi)*sinc(wc*(k-0.5*M)/pi);h=hd＇.*w; omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=abs(mag); plot(omega/pi,magdb,＇b＇);grid;w=ws-wp; M1=32;k2=-M1:M1;wc=(wp+ws)/2; hd=sinc(w*k2/2).*(sin(wc*k2)./(k2.*pi));hd(M1+1)=wc/pi; omega2=linspace(0,pi,512);mag2=freqz(hd,[1],omega2);magdb2=abs(mag2); hold on; plot(omega2/pi,magdb2,＇r＇); legend(＇逐步衰减＇,＇过渡带渐变＇); title(＇M1=32＇) 【研讨题目】 中等题 4．频率取样法FIR 数字滤波器 （1）（M5-6）利用频率取样法设计某I型线性相位带通FIR滤波器，其通带截频分别为 Wp1=0.3p rad, Wp2=0.5p rad （2）（M5-7）在通带和阻带间增加1个过渡点，重新设计该滤波器。过渡点的最佳幅度由实验确定。

【设计步骤】 （1） 学习利用频率取样法设计某I型线性相位带通FIR滤波器的方法；

（2） 确定理想滤波器的幅度函数；

（3） 确定理想滤波器的相位；

（4） 用实验法确定满足要求的滤波器的阶数；

（5） 用实验法确定过度点。

【仿真结果】 【结果分析】 滤波器的阻带衰减，滤波器的阶数与设计结果的关系。

滤波器的阶数越高，滤波器幅度响应越接近理想滤波器，滤波效果更好，不影响阻带衰减。滤波器阶数高也就使计算的难度增大，在能够满足工程需要的前提下阶数越小越好。

【自主学习内容】 【阅读文献】 【发现问题】 (专题研讨或相关知识点学习中发现的问题)：
在通带和阻带间增加过渡点的值比较难确定，找一个合适的方法比较重要。

【问题探究】 【问题探究】 如何确定过渡点的最佳幅度。

在确定最佳幅度时我采用的是类似于程序设计中的二分检索的方法，先设两个加大的值，以倍数的方式逐次向理想值靠近。这样可以加快速度的找出理想值。

【仿真程序】 (1) M=64;Wp1=0.3*pi;Wp2=0.5*pi;m=0:M/2; Wm=2*pi*m./(M+1); Ad=double([(Wm>=Wp1)&(Wm<=Wp2)]); Hd=Ad.*exp(-j*0.5*M*Wm); Hd=[Hd conj(fliplr(Hd(2:M/2+1)))]; h=real(ifft(Hd)); w=linspace(0.1,pi,1000); H=freqz(h,[1],w); plot(w/pi,20*log10(abs(H))); grid on; title(＇M=64带通滤波器＇); (2) M=64;Wp1=0.3*pi;Wp2=0.5*pi; Wp=(Wp1+Wp2)/2;m=0:M/2; Wm=2*pi*m./(M+1); mtr1=ceil(Wp1*(M+1)/(2*pi)); mtr2=floor(Wp2*(M+1)/(2*pi))+2; Ad=double([(Wm>=Wp1)&(Wm<=Wp2)]); Ad(mtr1)=0.40; Ad(mtr2)=0.42; Hd=Ad.*exp(-j*0.5*M*Wm); Hd=[Hd conj(fliplr(Hd(2:M/2+1)))]; h=real(ifft(Hd)); w=linspace(0.1,pi,1000); H=freqz(h,[1],w); plot(w/pi,20*log10(abs(H))); grid on; title(＇过渡点T1=0.40，T2=0.42＇); 【研讨题目】 中等题 5．设计幅度响应逼近下图所示的数字高通滤波器，其中Wc=0.5p。要求：
(1) 将该数字高通滤波器设计成IIR数字滤波器，具体要求：分别设计成BW型、CBI型、CBII型和椭圆型滤波器，并比较设计结果。

(2) 将该数字高通滤波器设计成FIR数字滤波器，具体要求：
a. 采用窗函数法，分别利用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗截断，并将设计结果进行分析比较。

b. 采用频率取样法，讨论过渡点对滤波器阻带衰减的影响。

c. 采用Parks-McClellan算法。

d. 比较窗函数法、频率取样法和Parks-McClellan算法所设计的FIR滤波器。

(3) 所设计的IIR和FIR数字滤波器应具有基本相同的幅度响应。根据设计结果，从幅度响应、相位响应、滤波器阶数等对两类滤波器进行分析比较。

【数字高通滤波器设计指标】 , ，Ap=1dB，As=50dB 【设计步骤】 模拟低通滤波器 模拟低通滤波器设计指标 模拟滤波器高通设计指标 数字高通滤波器设计指标 (1) → → → → 模拟高通滤波器 数字高通滤波器 → （2） a 单位脉冲响应 滤波器的频率响应 数字滤波器的设计指标 确定过渡带、窗函数、函数类型 IDTFT h[k] 加窗 b 抽样后的频率响应 H[m] 滤波器的频率响应 数字滤波器的设计指标 确定过渡带、窗函数、函数类型 抽样 h[k] DTFT 【仿真结果】 (2)a. 矩形窗 汉纳窗 汉明窗 Blackman窗 凯泽窗 B． C Parks-McClellan算法 【结果分析】 （1） BW型需要的阶数最高为11，椭圆的阶数最少5，其中chey12型的阶数均为6.因而BW型稳定性好，成本高，虽然椭圆阶数低，但稳定性差，难以实现。

（2） a由图我们可以看出，从矩形窗、汉纳窗、哈明窗、布莱克曼窗的顺序，在相同截取长度下，衰减逐渐增加，同时相应的近似过渡带也逐渐加宽。可以看出在增加衰减的同时会增加过渡带的宽度。

b加频率抽样点数可以是阻带的衰减明显增加，而且随着抽样点数的增加，过渡带的宽度减小，衰减幅度增加 cParks-McClellan算法设计出的滤波器阻带的衰减是等波纹的 d．两种滤波器幅度响应差不多，但是FIR可设计成线性相位，IIR相位是非线性的。从阶数上来看可以看出IIR的阶数可以更小。

【自主学习内容】 凯泽窗，频率取样法，Parks-McClellan算法的matlab表示。

【阅读文献】 《数字信号处理》 【仿真程序】 （1）Ap=1;As=50;omegap=0.6*pi;omegas=0.4*pi;T=2;w0=1; wp=tan(omegap/2);ws=tan(omegas/2); wp1=w0/wp;ws1=w0/ws; [N,wc,]=buttord(wp1,ws1,Ap,As,＇s＇) [num,den]=butter(N,wc,＇s＇); [numt,dent]=lp2hp(num,den,w0); [numd,dend]=bilinear(numt,dent,0.5); w=linspace(0,pi,512); h=freqz(numd,dend,w); norm=max(abs(h)); numd=numd/norm; subplot(121);plot(w/pi,abs(h)); subplot(122);plot(w/pi,20*log10(abs(h)/norm)); w=[omegap omegas]; h=freqz(numd,dend,w); A=20*log10(abs(h)/norm) 其他三个类似，不再赘述。

（2） omegap=0.6*pi;omegas=0.4*pi;Ap=1;As=50; N=ceil(1.8*pi/(omegap-omegas)); N=mod(N+1,2)+N M=N-1;w=ones(1,N);omegac=0.5*pi;k=0:M; hd=-(omegac/pi)*sinc(omegac*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1; h=hd.*w; omega=linspace(0,pi,512); mag=freqz(h,[1],omega); plot(omega/pi,20*log10(abs(mag)));grid;xlabel(＇Normalized frequency＇);ylabel(＇Gain in dB＇);（其他三个类似） 凯泽窗 Ap=1;As=50;Rp=1-10.^(-0.05*Ap);Rs=10.^(-0.05*As);f=[0.4,0.6];a=[0,1];dev=[Rp,Rs]; [M,Wc,beta,ftype] = kaiserord(f,a,dev); M=mod(M,2)+M,h = fir1(M,Wc,ftype,kaiser(M+1,beta)); omega=linspace(0,pi,512);mag=freqz(h,[1],omega); plot(omega/pi,20*log10(abs(mag))); grid; xlabel(＇Normalized frequency＇);ylabel(＇Gain,dB＇); （3）Fp=0.6;Fs=0.4;ds=0.0017;dp=ds; f=[Fs Fp];a=[0 1];dev=[ds dp]; [M,fo,ao,w] = remezord(f,a,dev); h = remez(M,fo,ao,w); w=linspace(0,pi,512); mag=freqz(h,[1],w); plot(w/pi,20*log10(abs(mag))); xlabel(＇Normalized frequency＇); ylabel(＇Gain, dB＇);grid;