DSP研究性学习报告滤波器设计.docx
《DSP研究性学习报告滤波器设计.docx》由会员分享，可在线阅读，更多相关《DSP研究性学习报告滤波器设计.docx（18页珍藏版）》请在冰豆网上搜索。
DSP研究性学习报告滤波器设计
《数字信号处理》课程研究性学习报告
数字滤波器设计专题研讨
【目的】
（1）掌握IIR和FIR数字滤波器的设计方法及各自的特点。
（2）掌握各种窗函数的时频特性及对滤波器设计的影响。
（3）培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。
【研讨题目1】
设计一个IIR数字低通滤波器，其能取代下列指标的模拟低通滤波器（系统的抽样频率为44.1kHz）
fp=2kHz，fs=10kHz,Ap=0.5dB,As=50dB
（1）分别用双线性变换和冲激响应不变法设计一个BW型数字低通滤波器，并进行比较。
（2）用双线性变换分别设计ChebyshevI型ChebyshevII型和椭圆型数字低通滤波器，并进行比较。
【温磬提示】
在数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。
但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。
【设计步骤】
由已知数字滤波器指标确定相应的模拟滤波器指标，然后设计原型模拟滤波器，再将其变换为相应的数字滤波器
即：
【结果分析】
【仿真结果】
双线性不变法
Ap=0.2445
As=50.0000
脉冲响应不变法：
Ap=0.4618
As=60.4050
ChebyshevI型
Ap=0.4999
As=71.0563
ChebyshevII型
Ap=0.5000
As=53.0577
椭圆型滤波器
Ap=0.5000
As=51.2355
【结果分析】
双线性变换和冲激响应不变法所设计的滤波器的性能有什么不同。
BW型、ChebyshevI型、ChebyshevII型和椭圆型滤波器的零极点分布各有什么特点
由图得脉冲响应不变法的主要优点是模拟频率与数字频率之间的关系是线性的，主要缺点是存在频谱混叠，使得阻带衰减不满足条件，不适合设计高通和带阻滤波器。
双线性变换法主要优点是避免了频谱混叠，其缺点是模拟频率与数字频率之间的关系是非线性的。
由Chebyshev I型、Chebyshev II型和椭圆型滤波器的仿真结果可以知道零极点分布的特点。
在相同的设计指标下，BW型滤波器的阶数最高，椭圆滤波器的阶数最低。
而阶数相等，它们的裕量也不同。
零极点分布都是关于虚轴对称且都在单位圆内。
总体来看，BW型滤波器最容易实现，而椭圆滤波器不易实现。
【自主学习内容】
利用matlab实现脉冲响应不变法和双线性变换法的程序语言
【阅读文献】
《数字信号处理》陈后金
【问题探究】
研讨探究的是双线性变换法和脉冲响应不变法以及BW型、Chebyshev I型、
Chebyshev II型和椭圆型滤波器的特点；
双线性变换法的优缺点
优点：
不会产生频谱混叠
缺点：
数字滤波器和模拟滤波器的频率关系非线性
脉冲响应不变法的优缺点
优点：
数字滤波器和模拟滤波器的频率关系为线性
缺点：
存在频谱混叠，不能用脉冲响应不变法设计高通、带阻等滤波器。
BW型、Chebyshev I型、Chebyshev I I型和椭圆型滤波器的特点;
1. 在相同的设计指标下，一般来说，BW型滤波器的阶数最高，椭圆滤波器的阶数最低。
即使阶数相等，它们的裕量也不同。
2．在滤波器的实现过程中，
BW型滤波器最容易实现，而椭圆滤波器不易实现（因为它的系统函数H（s）的极点离
jw轴最近）。
Cheby2型做的幅度响应在ws之后没有波动
【仿真程序】
双线性变换法
>>clear;
>>FS=44100;fp=2000;fs=10000;
>>Ap=0.5;As=50;
>>wp=fp*2*pi;ws=fs*2*pi;wp1=wp/FS;ws1=ws/FS;
>>OmegaP=2*FS*tan（wp1/2）;
>>OmegaS=2*FS*tan（ws1/2）;
>>[N,wc]=buttord（OmegaP,OmegaS,Ap,As,'s'）;
>>[bt,at]=butter（N,wc,'s'）;
>>[bz,az]=bilinear（bt,at,FS）;
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=1.167808e-024.
>Inbilinearat89
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=1.167808e-024.
>Inbilinearat90
>>w=linspace（0,pi,512）;h=freqz（bz,az,w）;norm=max（abs（h））;
>>bz=bz/norm;
>>plot（（w/pi）,20*log10（abs（h）/norm））;
>>xlabel（'Ω／π'）;ylabel（'幅值'）;
>>w=[wp1ws1];h=freqz（bz,az,w）;
>>fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
Ap=0.2445
>>fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
As=50.0000
>>gridon;
脉冲响应不变法
>>clear;
>>FS=44100;
>>fp=2000;fs=10000;
>>ws=fs*2*pi;
>>wp=fp*2*pi;
>>Ws=ws/FS;Wp=wp/FS;
>>Ap=0.5;As=50;
>>N=buttord（wp,ws,Ap,As,'s'）;
>>fprintf（'N=%.0f\n',N）;
N=10
>>wc=wp/10^（0.1*Ap-1）^（1/N/2）;
>>[n,d]=butter（N,wc,'s'）;
correct/robust.
CoeffsofB（s）/A（s）arereal,butB（z）/A（z）hascomplexcoeffs.
Probablecauseisrootingofhigh-orderrepeatedpolesinA（s）.
>Inimpinvarat122
>>w=linspace（0,pi,512）;
>>h=freqz（n,d,w）;
>>norm=max（abs（h））;
>>n=n/norm;
>>plot（w/pi,20*log10（abs（h）/norm））;
>>xlabel（'Ω／π'）;ylabel（'幅值'）;
>>w=[WpWs];
>>h=freqz（n,d,w）;
>>fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
Ap=0.4618
As=60.4050
>>gridon;
Chebyshev I型
>>clear;
>>FS=44100;
>>fp=20000;fs=10000;
>>Ap=0.5;As=50;
>>wp=fp*2*pi;
>>ws=fs*2*pi;
>>wp1=wp/FS;
>>ws1=ws/FS;
>>OmegaP=2*FS*tan（wp1/2）;
>>OmegaP=2*FS*tan（ws1/2）;
>>[N,wc]=cheb1ord（OmegaP,OmegaS,Ap,As,'s'）;
>>[bt,at]=cheby1（N,Ap,wc,'s'）;
>>[bz,az]=bilinear（bt,at,FS）;
>>w=linspace（0,pi,512）;
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=9.975559e-023.
>Inbilinearat89
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=9.975559e-023.
>Inbilinearat90z
>>h=freqz（bz,az,w）;
>>norm=max（abs（h））;
>>bz=bz/norm;subplot（2,1,1）;
>>plot（（w/pi）,20*log10（abs（h）/norm））;
>>xlabel（'Ω／π'）;ylabel（'幅值'）;
>>[r,p,k]=residuez（bz,az）;
>>subplot（2,1,2）;
>>zplane（bz,az）;
>>w=[wp1ws1];
>>h=freqz（bz,az,w）;
>>fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;
Ap=0.4999
>>fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
As=71.0563
Chebyshev I I型
>>clear;
>>FS=44100;
>>fp=20000;fs=10000;
>>Ap=0.5;As=50;
>>wp=fp*2*pi;
>>ws=fs*2*pi;
>>wp1=wp/FS;
>>ws1=ws/FS;
>>OmegaP=2*FS*tan（wp1/2）;
>>OmegaP=2*FS*tan（ws1/2）;
>>[N,wc]=cheb2ord（OmegaP,OmegaS,Ap,As,'s'）;
>>[bt,at]=cheby2（N,As,wc,'s'）;
>>[bz,az]=bilinear（bt,at,FS）;
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=1.956853e-024.
>Inbilinearat89
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=1.956853e-024.
>Inbilinearat90
>>w=linspace（0,pi,512）;
>>h=freqz（bz,az,w）;
>>norm=max（abs（h））;
>>bz=bz/norm;
>>subplot（2,1,1）;
>>plot（（w/pi）,20*log10（abs（h）/norm））;
>>xlabel（'Ω／π'）;ylabel（'幅值'）;
>>[r,p,k]=residuez（bz,az）;
>>subplot（2,1,2）;
>>zplane（bz,az）;
>>w=[wp1ws1];
>>h=freqz（bz,az,w）;
>>fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
Ap=0.5000
As=53.0577
椭圆型滤波器
>>clear;
>>FS=44100;
>>fp=20000;fs=10000;
>>Ap=0.5;As=50;
>>wp=fp*2*pi;
>>ws=fs*2*pi;
>>wp1=wp/FS;
>>ws1=ws/FS;
>>OmegaP=2*FS*tan（wp1/2）;
>>OmegaP=2*FS*tan（ws1/2）;
>>[N,wc]=ellipord（OmegaP,OmegaS,Ap,As,'s'）;
>>[bt,at]=ellip（N,Ap,As,wc,'s'）;
>>[bz,az]=bilinear（bt,at,FS）;
>>w=linspace（0,pi,512）;
>>h=freqz（bz,az,w）;
>>norm=max（abs（h））;
>>bz=bz/norm;
>>subplot（2,1,1）;
>>plot（（w/pi）,20*log10（abs（h）/norm））;
>>xlabel（'Ω／π'）;ylabel（'幅值'）;
>>[r,p,k]=residuez（bz,az）;
>>subplot（2,1,2）;zplane（bz,az）;
>>w=[wp1ws1];
>>h=freqz（bz,az,w）;
>>fprintf（'Ap=%.4f\n',-20*log10（abs（h
（1））））;>>fprintf（'As=%.4f\n',-20*log10（abs（h
（2））））;
Ap=0.5000
As=51.2355
【研讨题目2】
用窗函数法设计FIR数字高通滤波器，分别利用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗截断。
讨论用窗函数法设计FIR数字高通滤波器时如何确定滤波器的指标，比较相同过渡带时用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗设计滤波器的阶数。
【温馨提示】
窗函数法设计FIR滤波器的基本思想是在时域逼近理想滤波器的单位脉冲响应。
首先根据待逼近的理想滤波器的频率响应Hd（ej），由IDTFT求出理想滤波器的单位脉冲响应hd[k]，再将无限长的hd[k]加窗截断得到有限长序列h[k]。
为了获得线性相位FIR滤波器，在窗函数法设计FIR滤波器的过程中，需要将线性相位因子ej（0.5M加入理想滤波器的频率响应Hd（ej）。
常用窗函数除矩形窗外，还有Hann（汉纳）窗、Hamming（哈明）窗、Blackman（布莱克曼）窗、Kaiser（凯泽）窗等，这些窗函数的在频域的特征（主瓣宽度、旁瓣幅度）不同，使得所设计的滤波器过渡带宽度和阻带衰减也不同。
【设计步骤】
设计了下列指标的线性相位FIR高通滤波器：
Wp=0.67π,Ws=0.53π,Ap=0.3dB,As=50dB
使用不同的窗进行加窗截断：
采用矩形窗加窗截断
采用Hamming窗加窗截断，
采用blackman窗加窗截断
采用hanning窗加窗截断
采用凯泽窗加窗截断
【仿真结果】
N=51
采用矩形窗加窗截断
Ap≈0dBAs≈21dB
采用Hamming窗加窗截断
Ap≈0dBAs≈54 dB
采用blackman窗加窗截断
Ap≈0dBAs≈75dB
采用hanning窗加窗截断
Ap≈0dBAs≈44dB
采用凯泽窗加窗截断
Ap≈0dBAs≈52dB
【结果分析】
矩形窗设计出的FIR滤波器阻带衰减最小,布莱克曼窗阻带衰减最大
矩形窗、汉纳窗、哈明窗、布莱克曼窗比较
窗的类型
主瓣宽度
近似过度宽度
Ap
As
矩形
0.09
0.82
21
Hann
0.0064
0.056
44
Hamming
0.0022
0.019
53
Blackman
0.0002
0.0017
74
【自主学习内容】
用matlab实现设计窗函数，部分窗函数的调用形式为
W=hanning（N）
W=hamming（N）
W=blackman（N）
W=baiser（N）
其中N是窗函数的参数，返回的w是一个长度为N的列向量，给出了窗函数N点的取值。
【阅读文献】
《数字信号处理》陈后金
【问题探究】
矩形窗的过渡带最窄，但设计出的FIR滤波器的阻带衰减最小。
利用布莱克曼窗设计出的FIR滤波器阻带衰减最大，但其过渡带也最宽。
所以过渡带宽度与阻带衰减是一对矛盾。
工程应用中，在满足阻带衰减的前提下，尽可能地选择主瓣宽度较小的窗函数，对二者进行折中权衡处理。
【仿真程序】
矩形窗
>>Wp=0.67*pi;Ws=0.53*pi;Ap=0.3;As=50;
>>N=ceil（7*pi/（Wp-Ws））;
>>N=mod（N+1,2）+N;
>>M=N-1;
>>fprintf（'N=%.0f\n',N）;
N=51
>>w=boxcar（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）;magdb=20*log10（abs（mag））;
>>plot（omega/pi,magdb）;
哈明窗：
>>Wp=0.67*pi;Ws=0.53*pi;Ap=0.3;As=50;
>>N=ceil（7*pi/（Wp-Ws））;
>>N=mod（N+1,2）+N;
>>M=N-1;fprintf（'N=%.0f\n',N）;
N=51
>>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）;magdb=20*log10（abs（mag））;plot（omega/pi,magdb）
布莱克曼窗
>>Wp=0.67*pi;Ws=0.53*pi;Ap=0.3;As=50;
>>N=ceil（7*pi/（Wp-Ws））;N=mod（N+1,2）+N;
>>M=N-1;
>>fprintf（'N=%.0f\n',N）;
N=51
>>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;
>>hd（0.5*M+1）=hd（0.5*M+1）+1;
>>mag=freqz（h,[1],omega）;magdb=20*log10（abs（mag））;
>>plot（omega/pi,magdb）;
hanning窗
>>Wp=0.67*pi;Ws=0.53*pi;Ap=0.3;As=50;
>>N=ceil（7*pi/（Wp-Ws））;N=mod（N+1,2）+N;
>>
>>M=N-1;
>>fprintf（'N=%.0f\n',N）;
N=51
>>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）;magdb=20*log10（abs（mag））;
>>plot（omega/pi,magdb）;
凯泽窗
>>Ap=0.3;As=50;
>>Rp=1-10.^（-0.05*Ap）;Rs=10.^（-0.05*As）;
>>f=[0.53,0.67];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）））;
升级会员 