03DSP研究性学习报告数字滤波器设计.docx
《03DSP研究性学习报告数字滤波器设计.docx》由会员分享，可在线阅读，更多相关《03DSP研究性学习报告数字滤波器设计.docx（35页珍藏版）》请在冰豆网上搜索。
03DSP研究性学习报告数字滤波器设计
《数字信号处理》课程研究性学习报告
姓名
学号
同组成员
指导教师
时间
数字滤波器设计
【研讨题目】基本题
1．IIR数字滤波器设计
设计一个IIR数字低通滤波器，其能取代下列指标的模拟低通滤波器（系统的抽样频率为44.1kHz）
fp=2kHz，fs=10kHz,Ap=0.5dB,As=50dB
（1）分别用双线性变换和冲激响应不变法设计一个BW型数字低通滤波器，并进行比较。
（2）用双线性变换分别设计ChebyshevI型ChebyshevII型和椭圆型数字低通滤波器，并进行比较。
【温磬提示】
在数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。
但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。
【设计步骤】
一、脉冲响应不变法
1.将数字滤波器的频率指标{
}转换为模拟滤波器的频率指标{wk}
2.由模拟滤波器的指标设计模拟滤波器的H（s）。
模拟带通滤波器的设计步骤：
（1）由带通滤波器的上下截频确定变换式中的参数
（2）确定原型低通滤波器的通带截频
（3）设计通带截频为1（rad/s）、阻带截频为
、通带衰减为ApdB、阻带衰减为AsdB的原型低通滤波器
（4）将原型低通滤波器转换为带通滤波器HBP（s）
3.利用脉冲响应不变法，将H（s）转换H（z）。
二、双线性变换法
1.将数字滤波器的频率指标{
}转换为模拟滤波器的频率指标{wk}
2.由模拟滤波器的指标设计模拟滤波器的H（s）。
模拟带通滤波器的设计步骤：
（1）由带通滤波器的上下截频确定变换式中的参数
（2）确定原型低通滤波器的通带截频
（3）设计通带截频为1（rad/s）、阻带截频为
、通带衰减为ApdB、阻带衰减为AsdB的原型低通滤波器
（4）将原型低通滤波器转换为带通滤波器HBP（s）
3.利用双线性变换法，将H（s）转换H（z）。
【仿真结果】
所设计滤波器的幅度响应和相位响应
BW型、ChebyshevI型、ChebyshevII型和椭圆型滤波器的零极点分布
【结果分析】
双线性变换和冲激响应不变法所设计的滤波器的性能有什么不同。
BW型、ChebyshevI型、ChebyshevII型和椭圆型滤波器的零极点分布各有什么特点。
脉冲响应不变法：
clearall;
wp=2000*2*pi;
ws=10000*2*pi;
Fs=44100;
Wp=wp/Fs;Ws=ws/Fs;Ap=0.5;As=50;
N=buttord（wp,ws,Ap,As,'s'）;
fprintf（'N=%.0f\n',N）
wc=wp/（10^（0.1*Ap）-1）^（1/N/2）;
[numa,dena]=butter（N,wc,'s'）;
[numd,dend]=impinvar（numa,dena,Fs）;
w=linspace（0,pi,2048）;
h=freqz（numd,dend,w）;
norm=max（abs（h））;
numd=numd/norm;
figure
（1）;zplane（numa,dena）;
figure
（2）;plot（w/pi,20*log10（abs（h/norm）））;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
w=[WpWs];
h=freqz（numd,dend,w）;
fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
gridon;
双线性变化法：
clearall;
wp=2000*2*pi;ws=10000*2*pi;
Fs=44100;
Wp=wp/Fs;Ws=ws/Fs;Ap=0.5;As=50;
N=buttord（wp,ws,Ap,As,'s'）;
fprintf（'N=%.0f\n',N）
wc=wp/（10^（0.1*Ap）-1）^（1/N/2）;
[numa,dena]=butter（N,wc,'s'）;
[numd,dend]=bilinear（numa,dena,Fs）;
w=linspace（0,pi,2048）;
h=freqz（numd,dend,w）;
norm=max（abs（h））;
numd=numd/norm;
figure
（1）;zplane（numa,dena）;
figure
（2）;plot（w/pi,20*log10（abs（h/norm）））;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
w=[WpWs];
h=freqz（numd,dend,w）;
fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
gridon;
Chebyshev1：
clearall;
wp=2000*2*pi;ws=10000*2*pi;
Fs=44100;
Wp=wp/Fs;Ws=ws/Fs;Ap=0.5;As=50;
N=cheb1ord（wp,ws,Ap,As,'s'）;
fprintf（'N=%.0f\n',N）
wc=wp/（10^（0.1*Ap）-1）^（1/N/2）;
[numa,dena]=cheby1（N,Ap,wc,'s'）;
[numd,dend]=bilinear（numa,dena,Fs）;
w=linspace（0,pi,2048）;
h=freqz（numd,dend,w）;
norm=max（abs（h））;
numd=numd/norm;
figure
（1）;zplane（numa,dena）;
figure
（2）;plot（w/pi,20*log10（abs（h/norm）））;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
w=[WpWs];
h=freqz（numd,dend,w）;
fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
gridon;
Chebyshev2：
clearall;
wp=2000*2*pi;ws=10000*2*pi;
Fs=44100;
Wp=wp/Fs;Ws=ws/Fs;Ap=0.5;As=50;
N=cheb2ord（wp,ws,Ap,As,'s'）;
fprintf（'N=%.0f\n',N）
wc=wp/（10^（0.1*Ap）-1）^（1/N/2）;
[numa,dena]=cheby2（N,As,wc,'s'）;
[numd,dend]=bilinear（numa,dena,Fs）;
w=linspace（0,pi,2048）;
h=freqz（numd,dend,w）;
norm=max（abs（h））;
numd=numd/norm;
figure
（1）;zplane（numa,dena）;
figure
（2）;plot（w/pi,20*log10（abs（h/norm）））;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
w=[WpWs];
h=freqz（numd,dend,w）;
fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
gridon;
椭圆：
clearall;
wp=2000*2*pi;ws=10000*2*pi;
Fs=44100;
Wp=wp/Fs;Ws=ws/Fs;Ap=0.5;As=50;
N=ellipord（wp,ws,Ap,As,'s'）;
fprintf（'N=%.0f\n',N）
wc=wp/（10^（0.1*Ap）-1）^（1/N/2）;
[numa,dena]=ellip（N,Ap,As,wc,'s'）;
[numd,dend]=bilinear（numa,dena,Fs）;
w=linspace（0,pi,2048）;
h=freqz（numd,dend,w）;
norm=max（abs（h））;
numd=numd/norm;
figure
（1）;zplane（numa,dena）;
figure
（2）;plot（w/pi,20*log10（abs（h/norm）））;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
w=[WpWs];
h=freqz（numd,dend,w）;
fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
gridon;
【研讨题目】基本题
2．窗函数研究
分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性，并进行比较。
利用I型线性相位滤波器设计满足下列指标的FIR高通滤波器
Ωp=0.8π,Ωs=0.7π,Ap=0.3dBAs=40dB
【结果分析】
（1）Ap≈0dB
As≈20dB
（2）Ap=≈0.034dB
As≈44dB
（3）Ap≈0dB
As≈55dB
（4）Ap≈0dB
As≈78dB
（5）Ap≈0.057dB
As≈42dB
各种窗的比较
矩形窗
图为矩形窗设计的FIR滤波器在不连续点附近的幅度函数A（Ω）
Ap=-20lg（1-dp）≈0.82dB,
As=-20lg（ds）≈21dB
汉纳窗
图为汉纳窗设计的FIR滤波器在不连续点附近的幅度函数A（Ω）
Ap≈0.056dB,As≈44dB
哈明窗
图为哈明窗设计的FIR滤波器在不连续点附近的幅度函数A（Ω）
Ap≈0.019dB,As≈53dB
布莱克曼窗
图为布莱克曼窗设计的FIR滤波器在不连续点附近的幅度函数A（Ω）
Ap≈0.0017dB，As≈74dB
矩形窗、汉纳窗、哈明窗、布莱克曼窗比较
凯泽窗
I0（x）可用幂级数表示为
β是一可调参数，通过改变β的值可以调节窗函数的形状。
式中A=-20lg（min{δp，δs}）
滤波器阶数M或凯泽窗的长度则可由下式估计
【仿真程序】
（1）
Wp=0.8*pi;Ws=0.7*pi;Ap=0.3;As=40;
N=ceil（6.2*pi/（Wp-Ws））;
N=mod（N+1,2）+N;
M=N-1;
w=1;
Wc=（Wp+Ws）/2;
k=0:
M;
hd=-（Wc/pi）*sinc（Wc*（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（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
（2）
clearall£»
Wp=0.8*pi;Ws=0.7*pi;Ap=0.3;As=40;
N=ceil（6.2*pi/（Wp-Ws））;
N=mod（N+1,2）+N;
M=N-1;
w=hanning（N）';
Wc=（Wp+Ws）/2;
k=0:
M;
hd=-（Wc/pi）*sinc（Wc*（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（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
（3）
Wp=0.8*pi;Ws=0.7*pi;Ap=0.3;As=40;
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）;
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（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
（4）
Wp=0.8*pi;Ws=0.7*pi;Ap=0.3;As=40;
N=ceil（11.4*pi/（Wp-Ws））;
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）;
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（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
（5）
Ap=0.3;As=40;
Rp=1-10.^（-0.05*Ap）;Rs=10.^（-0.05*As）;
f=[0.7,0.8];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（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
由结果分析可知，在矩形窗、汉纳窗、哈明窗、布莱克曼窗中，矩形窗的过渡带最窄，但利用它设计出的FIR滤波器的阻带衰减最小。
利用布莱克曼窗设计出的FIR滤波器阻带衰减最大，但其过渡带也最宽。
显然，减小了窗函数旁瓣的相对幅度却增加了其主瓣的宽度，即提高FIR滤波器阻带衰减是以增加过渡带宽度为代价的。
在工程应用中，在满足阻带衰减的前提下，尽可能地选择主瓣宽度较小的窗函数。
而在实际设计中，可由待设计的FIR数字滤波器阻带衰减或通带波动来确定窗函数的类型，有过渡带宽度估计窗函数的长度N（N=M+1）。
而凯泽窗是一种应用广泛的可调窗，它可以通过改变窗函数的形状来控制窗函数旁瓣的大小，而在设计中可根据滤波器的衰减指标来确定窗函数的形状。
【研讨题目】基本题
3．窗函数法设计FIR数字滤波器
（1）分别用Blackman窗和Kaiser窗法设计一个满足下列指标的线性相位的FIR低通滤波器
Ωp=0.4πrad,Ap=0.5dB,Ωs=0.6πrad,As=55dB
【设计步骤】
（1）FIR低通滤波器设计步骤如下：
a.根据所需设计的滤波器，确定线性相位滤波器的类型（Ⅰ型，Ⅱ型，Ⅲ型，Ⅳ型）；
b.确定理想滤波器的幅度函数
；
c.确定理想滤波器的相位
。
，对Ⅰ型和Ⅱ型线性相位FIR滤波器β=π；
d.由式
，计算
；
e.对
进行加窗截断。
这道题我们设计的是一个线性相位的FIR低通滤波器
1,首先我们设计用Blackman窗设计一个满足下列指标的线性相位的FIR低通滤波器
选取理想低通滤波器的截频为
用Blackman窗可使滤波器满足指标。
由过渡带宽度得滤波器长度N需满足
用Ⅰ型滤波器，取N=53，由式
得
用长度为N=53的Blackman窗截断
即得所要求的FIR低通滤波器的
。
2,接下来，我们用Kaiser窗函数法设计一个满足下列指标的线性相位的FIR低通滤波器
（1）由给定指标确定待逼近理想低通的截频Wc
由于理想低通滤波器的|H（ejW）|在截频Wc处收敛于0.6，因此常将截频Wc取在过渡带的中点
（2）由给定指标确定Kaiser窗的参数N和b
所以A=-20lg（min{
}）=As=45dB
又因为I型线性相位滤波器阶数必须是偶数，取M=26，
所以
（3）设计截频Wc=0.5π的I型线性相位FIR低通滤波器
【仿真结果】
blackman图像
Kaiser图像
【结果分析】
比较两种窗的设计结果
（1）用Blackman窗设计的FIR低通滤波器N=53，通带和阻带衰减分别为Ap≈0dB，As≈78dB。
（2）用Kaiser窗函数法设计的线性相位FIR数字滤波器长度N=27，Kaiser窗的参数β=3.9754.滤波器通带和阻带衰减分别为Ap≈0dB，As≈46dB
【问题探究】
通过实验讨论如何控制滤波器的阻带衰减
比较分别用Blackman窗和Kaiser窗设计的两种结果可知，所设计出的滤波器都满足设计指标。
相比于常用窗函数，用Kaiser窗设计出的滤波器阶数较低，但滤波器的阻带波纹衰减较慢。
【仿真程序】
（1）Blackman窗
Wp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=55;
N=ceil（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）;
plot（omega/pi,20*log10（abs（mag）））;
grid;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
（2）Kaiser窗
wp=0.4*pi;ws=0.6*pi;As=55;Ap=0.5;
M=ceil（（As-7.95）/（ws-wp）/2.285）
M=M+mod（M,2）
beta=0.1102*（As-8.7）;
w=kaiser（M+1,beta）;
wc=（wp+ws）/2;
alpha=M/2;
k=0:
M;
hd=（wc/pi）*sinc（（wc/pi）*（k-alpha））;h=hd.*w';
omega=linspace（0,pi,512）;
mag=freqz（h,[1],omega）;
magdb=20*log10（abs（mag））;
plot（omega/pi,magdb）;
grid;
xlabel（'Normalizedfrequency'）;
ylabel（'Gain,dB'）;
axis（[0,1,-70,0]）;
（2）（M5-5）在用窗口法设计FIR滤波器时，由于理想滤波器的幅度响应在截频处发生突变，使得设计出的滤波器的幅度响应发生振荡，这个现象被称为Gibbs现象。
解决这个问题的一个方案是本书中介绍的用逐步衰减的窗函数。
另一个方案是使理想滤波器过渡带为渐变的，如下图所示具有线性过渡带的理想低通滤波器的频率响应，试用窗口法设计逼近该频率响应的FIR滤波器。
题3图
【
（2）单位脉冲响应证明】
试证该滤波器的单位脉冲响应为
其中：
，
【设计步骤】
首先用逐步衰减的窗函数（第一方案）设计一个FIR滤波器，再设计一个FIR滤波器，使其理想滤波器过渡带为渐变的，并用矩形窗截断（第二方案）。
然后分析两种方法设计出来的滤波器，得出结论。
【仿真结果】
渐变的窗函数选择hamming窗。
为了简便研究过程，设Ωp=0.55π、Ωs=0.45π、
=25dB、
=1dB。
设hamming窗阶数为M，矩形窗的长度为M1，图中蓝线为第一种方案涉及到滤波器，红线为第二种方案设计的滤波器。
易知在本题中M=7。
M1=5
M1=8时
M1=15时
M1=30
M1=80
【结果分析】
由仿真结果可知，第一种方案的过渡带明显短于第二种方案；当M=7时，第一种方案几乎可以完全消除Gibbs现象，但M1=8时，第二种方案仍然可以看到明显的通阻带波动，故消除Gibbs现象方案二需要的阶数更高。
【仿真程序】
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）;
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）;
holdon;
plot（omega2/pi,magdb2,'r'）;
【研讨题目】中等题
4．频率取样法FIR数字滤波器
（1）（M5-6）利用频率取样法设计某I型线性相位带通FIR滤波器，其通带截频分别为
Ωp1=0.3πrad,Ωp2=0.5πrad
（2）（M5-7）在通带和阻带间增加1个过渡点，重新设计该滤波器。
过渡点的最佳幅度由实验确定。
【设计步骤】
频率取样法FIR数字滤波
升级会员 