滤波器设计—fdatool工具及代码设计

1、 黎美琪 201300800610 通信工程2班实验一 基于FDATool的FIR及IIR滤波器设计1、 实验目的： 1.学会使用fdatool设计滤波器 2.分析比较不同滤波器的特性2、 实验条件 Matlab 2013b pc机三、实验过程 知识储备：使用FDATool设计FIR滤波器的具体步骤1.1 滤波器指标 若需要设计一个16阶的FIR滤波器（h(0)=0），给定的参数如下： (1) 低通滤波器 注：阶数问题 这个选项直接影响滤波器的性能，阶数越高，性能越好，但是相应在FPGA实现耗用的资源需要增多。在这个设置中提供2个选项：Specify order和Minimum order，S

2、pecify order是工程师自己确定滤波器的阶数，Minimum order是让工具自动确定达到期望的频率相应所需要的最小阶数。(2) 采样频率Fs为8kHz，滤波器Fpass为3kHz，Fstop为3050hz (3) 输入序列位宽为9位（最高位为符号位） 在此利用MATLAB来完成FIR滤波器系数的确定。 1.2 打开MATLAB的FDATool MATLAB集成了一套功能强大的滤波器设计工具FDATool（Filter Design & Analysis Tool），可以完成多种滤波器的设计、分析和性能评估。在MATLAB主命令窗口内键入“fdatool”，打开FDATool

3、程序界面，如图B.2所示。 1.3 滤波器设计 FDATool界面左下侧排列了一组工具按钮，其功能分别如下所述： 滤波器转换（TransForm Filter） 设置量化参数（Set Quantization Parameters） 实现模型（Realize Model） 导入滤波器（Import Filter） 多速率滤波器（Multirate Filter） 零极点编辑器（Pole-zero Editor） 设计滤波器（Design Filter） 选择其中的选择Design Filter按钮，进入设计滤波器界面，进行下列选择，如图B.3所示。 滤波器类型（Filer Type）为低通（L

4、ow Pass） 设计方法（Design Method）为FIR/IIR，分别采用Equiripple、Least-squares、Window、Constr.LeastPth-norm 、Constrained Equiripple、Constr.Band Equiripple（FIR滤波器设计）和 Butterworth、Chebyshev Type I、Chebyshev Type I、Elliptic、Maximally flat、Least Pth-norm、Constr.LeastPth-norm（IIR滤波器设计）。 滤波器阶数（Filter order）定制为15 Fs =8k

5、Hz，Fpass=3kHz ，Fstop=3050hz最后单击Design Filter 图标，让MATLAB 计算FIR 滤波器系数并作相关分析。 其系统函数H(z)可用下式来表示： 显然上式可以写成: 即可以看成是一个15 阶的FIR 滤波器的输出结果经过了一个单位延时单元，所以在FDATool 中,把它看成15 阶FIR 滤波器来计算参数。 1.4 滤波器分析 计算完FIR 滤波器系数以后，往往需要对设计好的FIR 滤波器进行相关的性能分析，以便了解该滤波器是否满足设计要求。分析操作步骤如下： 选择选择工具栏中“Magnitude and Phase Response”，启动幅、相频比较

6、分析；选择工具栏中的圆图图标，启动零极点位置分析（主要分析滤波器的稳定性，但是通过滤波器基本信息栏也可以直接看滤波器的稳定性）。除此之外，工具栏 的不同图标还可以对滤波器进行群延时、相位延时、单位脉冲响应、单位阶跃响应、滤波器的系数的分析（此实验中没有对各方面都展开分析）。在图形的左侧，显示了当下滤波器的基本信息，如：1.5 量化 可以看到，FDATool计算出的值是一个有符号的小数，如果建立的FIR滤波器模型需要一个整数作为滤波器系数，就必须进行量化，并对得到的系数进行归一化。为此，单击FDATool左下侧的工具按钮进行量化参数设置。量化参数有三种方式：双精度、单精度和定点。（注意：当使用D

7、SP处理器时，FTR滤波器的系数要量化，量化误差将导致所设计的滤波器形状发生很大的变化；对于IIR滤波器，会改变滤波器的特性，甚至导致滤波器不稳定。）1.6 导出滤波器系数 为导出设计好的滤波器系数，选择FDATool菜单的“File”“Export”命令，打开Export（导出）对话框，选择导出到工作区（Workplace）。这时FIR滤波器系数就存入到一个一维变量Num（名称按照顺序取为a-g）中了，不过这时Num中的元素是以小数形式出现的。对于IIR滤波器，导出得到矩阵SOS（名称按照顺序取为h、k-p）和G（名称按照顺序取为h1、k1-p1）。导出后的数据在workplace中的显示如

8、下：其中FIR滤波器系数滤波器名称 滤波器系数12345678Equiripple（a）-0.179-0.0570.086-0.1130.114-0.064-0.0760.587Least-squares（b）-0.0360.01310.027-0.0670.081-0.043-0.0870.59Window（c）-0.0370.01790.021-0.0640.083-0.048-0.0810.588Constr.LeastPth-norm （d）-0.0380.0170.022-0.0640.083-0.049-0.080.587Least Pth-norm（e）0.53420.2105-

9、0.48-0.255-0.07-0.259-0.102-0.06Constrained Equiripple（f）-0.0110.02710.024-0.0460.069-0.039-0.0830.588Constr.Band Equiripple（g）-0.179-0.0570.086-0.1130.114-0.064-0.0760.587滤波器名称 滤波器系数910111213141516Equiripple（a）0.5875-0.076-0.060.1142-0.110.0855-0.057-0.18Least-squares（b）0.5904-0.087-0.040.0814-0.07

10、0.02730.0131-0.04Window（c）0.588-0.081-0.050.0829-0.060.02140.0179-0.04Constr.LeastPth-norm （d）0.5872-0.08-0.050.0826-0.060.02250.017-0.04Least Pth-norm（e）-0.168-0.015-0.08-0.0890.008-0.1230.01650.307Constrained Equiripple（f）0.5881-0.083-0.040.0686-0.050.02450.0271-0.01Constr.Band Equiripple（g）0.5872

11、-0.076-0.060.1143-0.110.0855-0.057-0.18对于IIR滤波器：b,a=sos2tf(SOS);得到滤波器系统函数的分子分母多项式系数，只是差了一个比例因子k。K=cumprod(G); k=K(end);滤波输出结果filteredpWave可以这样求：filteredpWave=filter(b,a,pWave)*k; 导出的数与传输函数的系数的关系a，b含义，在MATLAB中函数如下公式：这种用法主要出现在信号处理方向，控制方向上主要使用z，而不是z-1，这是习惯问题，都一样。其次，了解了其具体函数，则它的传递函数就好写了，如下：h = tf(b,a,0.

12、1,'variable','z-1');     %0.1表示采样时间，你可以根据自己设计的滤波器来设定。h表示的传递函数就是上面公式中的H(z)。当然，你也可以用“Discrete Transfer Fcn”模块来建立传递函数，由于这个模块是直接以Z的形式，所以a,b系统应该倒过来写，相信这个你应该理解。 按照上述操作，得到的滤波器设计相关图形如下：FIR滤波器设计滤波器幅度、相位响应 滤波器零极点图形 滤波器单位脉冲响应Equiripple Least-squares Window（Kaiser Beta=0.5）滤波器幅度、相

13、位响应 滤波器零极点图形 滤波器单位脉冲响应Constr.LeastPth-normConstrained EquirippleConstr.Band EquirippleIIR滤波器设计滤波器幅度、相位响应 滤波器零极点图形 滤波器单位脉冲响应Butterworth（min order175）Chebyshev Type I（min order 32）Chebyshev Type II（min order 32）Elliptic（min order 12）Maximally flatLeastxPth-norm、Constr.Least Pth-norm虽然以上的IIR滤波器的极点都是在单位

14、圆上，但是上述滤波器都是稳定的，从其基本信息可以看出 图形数据分析： 1.通过观察FIR滤波器的“幅度相位图”及“脉冲响应图”可以发现，脉冲响应关于中点对称 的FIR滤波器无相位失真。对称的脉冲响应产生了通带内的线性相位关系。这样，线性相位关系也保证无相位失真；但是IIR滤波器的相位均是非线性的，即其不能保证无相位失真。 2.对于给定的滤波器形状，IIR滤波器较FIR滤波器所需要的系数要少。 3.从“零极点图”可以发现，滤波器的系数为N个，对应的有N-1个零点和N-1个极点，FIR滤波器的极点都位于单位圆内（原点），即所有的FIR滤波器都是稳定的。而IIR（它是递归滤波器）的极点由分母多项式所

15、确定，极点分布在单位圆内及单位圆上，这意味着不能保证其稳定性，但是从基本信息处可以看到上面所讨论的几种IIR滤波器是稳定的。 4.通常的几种模拟滤波器的特性：Butterworth滤波器在通带和阻带都是单调的，意味着它们在一个方向上平滑变化；Chebyshev Type I滤波器在阻带内是单调的，但是在通带内是有波纹的；Chebyshev Type II滤波器在通带内是单调的，但是在阻带内是有波纹的；Elliptic滤波器在通带和阻带内都有波纹。 5.对于同样的设计要求，切比雪夫滤波器(min order 32)比巴特沃斯滤波器(min order175)的阶数要低得多。补充说明：FIR和II

16、R比较小结FIR滤波器IIR滤波器H（n）有限长H（n）无限长极点固定在原点极点位于Z平面任意位置滤波器阶次高滤波器阶次低可严格的线性相位非线性相位一般采用非递归结构递归结构可以用FFT计算不能用FFT计算设计借助于计算机可用模拟滤波器设计可设计各种幅频特性和相频特性的滤波器可用于设计规格化的选频滤波器 实验二 利用函数法设计FIR及IIR滤波器1、 实验目的 1.熟悉掌握设计滤波器的操作步骤 2.学会用频率响应采样法和窗函数法设计FIR滤波器，并比较其性能 3.学会用单位脉冲响应法和双线性法设计IIR滤波器，并比较其性能2、 实验条件 Matlab 2013b pc机3、 实验过程滤波器指标

17、：Fs=8000hz；fpass=3000hz；fstop=3050；s=0.1；p=0.1 （1）窗函数法设计FIR低通滤波器（主要可以参见实验四）设计步骤1.在过渡带宽度的中间，选择通带边缘频率（Hz）： f1=fpass+(fstop-fpass)/2=30252.计算1=2*pi*f1/Fs,并将此值代入理想低通滤波器的脉冲响应h1n中： h1n=sin(n*1)/(n*pi)= sin(n*0.75625*pi)/(n*pi)3.计算滤波器阻带衰减为 Dbs=|20lg（0.1）|=20 通带边缘增益为 Dbp=20lg（1-0.1）=-0.91从表1中选择满足阻带衰减及其他滤波器要

18、求的窗函数矩形窗，用表中的公式计算所需窗的非零项目数目N=0.91*Fs/（fstop-fpass）=145.6选择奇数项145，这样脉冲响应完全可以对称，避免了滤波器产生相位失真，对于|n|<=(N-1)/2,设计窗函数wn。 表1 4.对于|n|<=(N-1)/2,从式hn=h1nwn计算（有限）脉冲响应，对于其他n值hn=0,此脉冲响应是非因果的。5.将脉冲响应右移(N-1)/2，确保第一个非零值在n=0处，使此低通滤波器为因果的。 编写代码如下：%LPF design 窗函数法clc;clear all;close all;omga_d=3025*pi/4000; %理想低

19、通滤波器的频率omga=0:pi/100:pi;N=145;w1=window(rectwin,N);%矩形窗 subplot(411);plot(1:N,w1);axis(1 N 0 1.1);%取值范围不同，矩形窗的中心位置不同title('窗函数时域形状');legend('rectwin');%标注窗函数名称M=floor(N/2);%M=72n=1:M;%hd=sin(n*omga_d)./(n*omga_d)*omga_d/pi;hd=sin(n*omga_d)./(n*pi);%理想低通滤波器对应的脉冲响应的右边序列hdd=fliplr(hd),1

20、/omga_d,hd; %绘制原脉冲函数图像，关于（N+1）/2对称subplot(412);plot(1:N,hdd);axis(1 N -0.2 0.45);title('脉冲函数hd（n）');h_d1=hdd.*w1'%被窗函数截取后的脉冲函数m=1:M;H_d1=2*cos(omga'*m)*h_d1(M+2:N)'+h_d1(M+1);%对截取后的脉冲函数做newdtft运算得到Hsubplot(413);plot(omga,H_d1);grid on;axis(0 pi 0 1);title('窗函数法得滤波器形状H');s

21、ubplot(414);plot(abs(fft(h_d1)');title('fft求得滤波器形状');图形显示：虽然矩形窗得到的低通滤波器有较好的特性，但是由于其边缘陡峭，旁瓣出现振铃，且变化较为明显。（2）频率响应采样法设计FIR低通滤波器。 设计原理及滤波器性能分析：.频率采样法是从频域出发，对给定的理想滤波器的频响 进行N点等间隔采样 然后以此Hd(k)作为实际FIR滤波器的频率特性采样值H(k)，即令： 由DFT定义，可以用这N个频域的采样值H(k)来惟一确定FIR的单位脉冲响应h(n)，即： 图示如下：对设计出的滤波器频率响应特性进行分析。由频域采样定理中

22、的内插公式可以知道，利用这N个频域采样值H(k)同样可以求得FIR滤波器的频率响应，这个将逼近理想滤波器的频响 。的内插公式为： （1）式中（w）是内插函数： （2）把（2）式带入（1）式，化简后可得： （3）从式(3)可以看到，在各频率采样点上，设计的滤波器，实际的频率响应严格地与理想滤波器的频率响应数值相等，即 但是在采样点之间的频率响应是由各采样点的加权内插函数叠加而形成的，因而有一定的逼近误差。该误差大小取决于理想频率响应的形状，理想频响特性变化越平缓，内插值越接近理想值，逼近误差越小；反之，如果采样点之间的理想频响特性变化越陡，则内插值与理想值之间的误差越大，因而在理想滤波器不连续点

23、的两边，就会产生尖峰，而在通带和阻带就会产生波纹。线性相位条件：FIR滤波器的最大优点是严格的线性相位特性。下面讨论为实现线性相位，在频域内采样得到的Hd(k)应满足什么条件。FIR滤波器具有线性相位的条件是h(n)是实序列，且满足h(n)=±h(N-1-n)，即h(n)关于对称，其中N为滤波器的长度。以第一类线性相位条件h(n)=h(N-1-n)(偶对称)为基础来推导频域采样Hd(k)满足的条件。理想滤波器的频率响应可以表示为： （4）为实现第一类线性相位条件，相位函数()和幅度函数Hg()应分别满足： （5） （6）在w=02*pi之间对进行N点等间隔采样，即，将w代入式（4）、

24、（5）、（6）并表示为k的函数，可以得到： （7） （8） （9）式(8)和式(9)就是频率采样值Hd(k)满足第一类线性相位的条件。设计步骤：1. 确定希望逼近的理想滤波器的频率响应2. 在频域内对进行k点等间隔采样，利用频率采样设计公式求频率采样值Hd(k)，采样间隔=k=1/512，这样在通带内共有c个采样点，分别是k=0，1，c-1。3. 用离散傅里叶逆变换求得要设计的实际滤波器的单位脉冲响应h(n)： 4. 根据傅里叶变换的定义求得实际滤波器的频率响应，验证是否满足滤波器技术指标的要求，主要验证滤波器的阻带衰减是否能够满足阻带的要求。Matlab代码中需要用到的函数知识>>

25、; help linspace linspace Linearly spaced vector. linspace(X1, X2) generates a row vector of 100 linearly equally spaced points between X1 and X2. linspace(X1, X2, N) generates N points between X1 and X2. For N = 1, linspace returns X2.编写代码如下：%LPF design 频率响应采样法%待设计滤波器的基本参数f_p=3000;%所设计滤波器的通带边缘频率f_s=

26、3050;%阻带边缘频率Fs=8000;%采样频率%All in Hzdelta_p=0.1;%通带波纹delta_s=0.1;%阻带波纹omg_p=2*pi*f_p/Fs;%通带边缘模拟频率转化为数字频率omg_s=2*pi*f_s/Fs;%阻带边缘模拟频率转化为数字频率K=512;%理想滤波器长度M=56;%单位脉冲响应单边长度C=floor(omg_p+omg_s)/2/pi*K);%理想滤波器中有值部分长度？Mag=ones(1,C),zeros(1,K-C);%幅度phi=linspace(0,pi,K)*M/2;%相位，范围为0-pi，间隔为pi/K，线性相位特性omg*M/2H=

27、Mag.*exp(-1j*phi);%|H|*ej*omg*M/2，待设计的理想滤波器H=H,conj(fliplr(H(1:K-1);%将H转换为因果的h0=fftshift(real(ifft(H);h=h0(K+M/2-M:K+M/2+M);%相移ej*omg*M/2,时域采样平移M/2个点subplot(211);plot(h);%脉冲响应axis(1 120 -0.2 0.8);title('脉冲函数hd（n）');subplot(212);plot(abs(fft(h);%所设计的滤波器的形状title('所设计的滤波器形状');图像显示：增加过渡带

28、抽样点，可加大阻带衰减，但是导致过渡带变宽增加N,使抽样点变密，减小过渡带带宽，但是增加了计算量优点：频域直接计算缺点：抽样频率只能是2/N或者/N的整数倍，截止频率不能任意取值（3）单位脉冲响应法设计IIR滤波器核心原理：通过对连续函数ha（t）等间隔采样得到离散序列ha（nT）。令h（n）=ha（nT），T为采样间隔。它是一种时域上的转换方法。设计步骤：1. 由模拟滤波器的传输函数Ha（s）求出其对应的模拟冲激响应ha（t）2. 对ha（t）每T秒进行一次采样，采样时间为t=nT，得到数字脉冲响应h（n）=ha（nT）3. 由h（n）进行z变换可以求出数字滤波器的传输函数H（z）4. 将z

29、用ej*Omga替换得到数字滤波器的频率响应H（Omga）编写代码如下：（用要求的参量得不到需要设计的滤波器，这里更换了参数）%LPF design 脉冲响应不变法(Butterworth)fp=100;fs=300;Rp=0.91;%通带衰减系数Rs=20;%阻带衰减系数Fs=1000;wp=2*pi*fp;ws=2*pi*fs;N,Wn=buttord(wp,ws,Rp,Rs,'s')%求滤波器所需阶数及截止频率B,A=butter(N,Wn,'s')%计算滤波器传输函数系数W=0:pi/100:pi;H,W=freqs(B,A,W);%滤波器频率响应%H=

30、20*log10(abs(H);plot(W,H),grid on;xlabel('rad');ylabel('幅度/db');title('所设计的滤波器形状'); axis(0 pi 0 1.1);运行结果N = 3 Wn = 876.3856B =1.0e+08 * 0 0 0 6.7311A = 1.0e+08 * 0.0000 0.0000 0.0154 6.7311图形显示：（4）双线性法设计IIR滤波器设计步骤1. 确定待求通带边缘频率fpass=3000hz，待求阻带边缘频率fstop=3050hz，待求阻带衰减20db，通带边缘

31、增益-0.91db。2. 用式子=2*f/Fs把由hz表示的待求边缘频率转化为由弧度表示的数字频率，得到pass、stop。 3. 计算预扭曲模拟频率以避免双线性变换带来的失真。由w=2Fs*tan（/2）求得wpass和wstop单位是弧度/秒。4. 计算所需滤波器的阶数5. 把wpass带入n阶模拟巴特沃斯滤波器传输函数H（s）中，并且对其进行双线性变换s=2Fs*（z-1）/（z+1）得到n阶数字传输函数H（z）。6. 由H(z)求出滤波器实现所需的差分方程。7. 把w=2Fs*tan（/2）代入即可得滤波器的形状|H（）|。编写代码如下：%LPF design 双线性变换法(Butte

32、rworth)%待设计滤波器的基本参数fpass=3000;fstop=3050;Fs=8000;%单位hzRs=0.1; %阻带波纹Rp=0.1; %通带波纹omg_pass=2*pi*fpass/Fs %通带截止频率数字频率omg_stop=2*pi*fstop/Fs %阻带截止频率数字频率，单位radwpass=2*Fs*tan(omg_pass/2)wstop=2*Fs*tan(omg_stop/2) %预扭曲模拟频率，单位rad/sn=ceil(log10(1/(Rs)2-1)/(2*log10(wstop/wpass) %所需滤波器阶数,ceil向上取整omg=0:pi/100:pi;H=1./sqrt(2*Fs*tan(omg./2)./wpass).(2*n)+1);plot(omg,abs(H);grid on;title('所设计的滤波器形状'); axis(0 pi 0