滤波器设计fdatool工具及代码设计详解.docx
《滤波器设计fdatool工具及代码设计详解.docx》由会员分享，可在线阅读，更多相关《滤波器设计fdatool工具及代码设计详解.docx（30页珍藏版）》请在冰点文库上搜索。
滤波器设计fdatool工具及代码设计详解
黎美琪201300800610通信工程2班
实验一基于FDATool的FIR及IIR滤波器设计
一、实验目的：
1.学会使用fdatool设计滤波器
2.分析比较不同滤波器的特性
二、实验条件
Matlab2013bpc机
三、实验过程
知识储备：
使用FDATool设计FIR滤波器的具体步骤
1.1滤波器指标
若需要设计一个16阶的FIR滤波器（h（0）=0），给定的参数如下：
（1）低通滤波器
注：
阶数问题
这个选项直接影响滤波器的性能，阶数越高，性能越好，但是相应在FPGA实现耗用的资源需要增多。
在这个设置中提供2个选项：
Specifyorder和Minimumorder，Specifyorder是工程师自己确定滤波器的阶数，Minimumorder是让工具自动确定达到期望的频率相应所需要的最小阶数。
（2）采样频率Fs为8kHz，滤波器Fpass为3kHz，Fstop为3050hz
（3）输入序列位宽为9位（最高位为符号位）在此利用MATLAB来完成FIR滤波器系数的确定。
1.2打开MATLAB的FDATool
MATLAB集成了一套功能强大的滤波器设计工具FDATool（FilterDesign&AnalysisTool），可以完成多种滤波器的设计、分析和性能评估。
在MATLAB主命令窗口内键入“fdatool”，打开FDATool程序界面，如图B.2所示。
1.3滤波器设计
FDATool界面左下侧排列了一组工具按钮，其功能分别如下所述：
滤波器转换（TransFormFilter）
设置量化参数（SetQuantizationParameters）
实现模型（RealizeModel）
导入滤波器（ImportFilter）
多速率滤波器（MultirateFilter）
零极点编辑器（Pole-zeroEditor）
设计滤波器（DesignFilter）
选择其中的选择DesignFilter按钮，进入设计滤波器界面，进行下列选择，如图B.3所示。
滤波器类型（FilerType）为低通（LowPass）
设计方法（DesignMethod）为FIR/IIR，分别采用Equiripple、Least-squares、Window、Constr.LeastPth-norm、ConstrainedEquiripple、Constr.BandEquiripple（FIR滤波器设计）和Butterworth、ChebyshevTypeI、ChebyshevTypeI、Elliptic、Maximallyflat、LeastPth-norm、Constr.LeastPth-norm（IIR滤波器设计）。
滤波器阶数（Filterorder）定制为15
Fs=8kHz，Fpass=3kHz，Fstop=3050hz
最后单击DesignFilter图标，让MATLAB计算FIR滤波器系数并作相关分析。
其系统函数H（z）可用下式来表示：
显然上式可以写成:
即可以看成是一个15阶的FIR滤波器的输出结果经过了一个单位延时单元，所以在FDATool中,把它看成15阶FIR滤波器来计算参数。
1.4滤波器分析
计算完FIR滤波器系数以后，往往需要对设计好的FIR滤波器进行相关的性能分析，以便了解该滤波器是否满足设计要求。
分析操作步骤如下：
选择选择工具栏中“MagnitudeandPhaseResponse”，启动幅、相频比较分析；选择工具栏中的圆图图标，启动零极点位置分析（主要分析滤波器的稳定性，但是通过滤波器基本信息栏也可以直接看滤波器的稳定性）。
除此之外，工具栏{
}的不同图标还可以对滤波器进行群延时、相位延时、单位脉冲响应、单位阶跃响应、滤波器的系数的分析（此实验中没有对各方面都展开分析）。
在图形的左侧，显示了当下滤波器的基本信息，如：
1.5量化
可以看到，FDATool计算出的值是一个有符号的小数，如果建立的FIR滤波器模型需要一个整数作为滤波器系数，就必须进行量化，并对得到的系数进行归一化。
为此，单击FDATool左下侧的工具按钮进行量化参数设置。
量化参数有三种方式：
双精度、单精度和定点。
（注意：
当使用DSP处理器时，FTR滤波器的系数要量化，量化误差将导致所设计的滤波器形状发生很大的变化；对于IIR滤波器，会改变滤波器的特性，甚至导致滤波器不稳定。
）
1.6导出滤波器系数
为导出设计好的滤波器系数，选择FDATool菜单的“File”→“Export”命令，打开Export（导出）对话框，选择导出到工作区（Workplace）。
这时FIR滤波器系数就存入到一个一维变量Num（名称按照顺序取为a-g）中了，不过这时Num中的元素是以小数形式出现的。
对于IIR滤波器，导出得到矩阵SOS（名称按照顺序取为h、k-p）和G（名称按照顺序取为h1、k1-p1）。
导出后的数据在workplace中的显示如下：
其中FIR滤波器系数
滤波器名称滤波器系数
1
2
3
4
5
6
7
8
Equiripple（a）
-0.179
-0.057
0.086
-0.113
0.114
-0.064
-0.076
0.587
Least-squares（b）
-0.036
0.0131
0.027
-0.067
0.081
-0.043
-0.087
0.59
Window（c）
-0.037
0.0179
0.021
-0.064
0.083
-0.048
-0.081
0.588
Constr.LeastPth-norm（d）
-0.038
0.017
0.022
-0.064
0.083
-0.049
-0.08
0.587
LeastPth-norm（e）
0.5342
0.2105
-0.48
-0.255
-0.07
-0.259
-0.102
-0.06
ConstrainedEquiripple（f）
-0.011
0.0271
0.024
-0.046
0.069
-0.039
-0.083
0.588
Constr.BandEquiripple（g）
-0.179
-0.057
0.086
-0.113
0.114
-0.064
-0.076
0.587
滤波器名称滤波器系数
9
10
11
12
13
14
15
16
Equiripple（a）
0.5875
-0.076
-0.06
0.1142
-0.11
0.0855
-0.057
-0.18
Least-squares（b）
0.5904
-0.087
-0.04
0.0814
-0.07
0.0273
0.0131
-0.04
Window（c）
0.588
-0.081
-0.05
0.0829
-0.06
0.0214
0.0179
-0.04
Constr.LeastPth-norm（d）
0.5872
-0.08
-0.05
0.0826
-0.06
0.0225
0.017
-0.04
LeastPth-norm（e）
-0.168
-0.015
-0.08
-0.089
0.008
-0.123
0.0165
0.307
ConstrainedEquiripple（f）
0.5881
-0.083
-0.04
0.0686
-0.05
0.0245
0.0271
-0.01
Constr.BandEquiripple（g）
0.5872
-0.076
-0.06
0.1143
-0.11
0.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.1,'variable','z^-1'）; %0.1表示采样时间，你可以根据自己设计的滤波器来设定。
h表示的传递函数就是上面公式中的H（z）。
当然，你也可以用“Discrete TransferFcn”模块来建立传递函数，由于这个模块是直接以Z的形式，所以a,b系统应该倒过来写，相信这个你应该理解。
按照上述操作，得到的滤波器设计相关图形如下：
FIR滤波器设计
滤波器幅度、相位响应滤波器零极点图形滤波器单位脉冲响应
Equiripple
Least-squares
Window（KaiserBeta=0.5）
滤波器幅度、相位响应滤波器零极点图形滤波器单位脉冲响应
Constr.LeastPth-norm
ConstrainedEquiripple
Constr.BandEquiripple
IIR滤波器设计
滤波器幅度、相位响应滤波器零极点图形滤波器单位脉冲响应
Butterworth（minorder175）
ChebyshevTypeI（minorder32）
ChebyshevTypeII（minorder32）
Elliptic（minorder12）
Maximallyflat
LeastxPth-norm、
Constr.LeastPth-norm
虽然以上的IIR滤波器的极点都是在单位圆上，但是上述滤波器都是稳定的，从其基本信息可以看出
图形数据分析：
1.通过观察FIR滤波器的“幅度相位图”及“脉冲响应图”可以发现，脉冲响应关于中点对称的FIR滤波器无相位失真。
对称的脉冲响应产生了通带内的线性相位关系。
这样，线性相位关系也保证无相位失真；但是IIR滤波器的相位均是非线性的，即其不能保证无相位失真。
2.对于给定的滤波器形状，IIR滤波器较FIR滤波器所需要的系数要少。
3.从“零极点图”可以发现，滤波器的系数为N个，对应的有N-1个零点和N-1个极点，FIR滤波器的极点都位于单位圆内（原点），即所有的FIR滤波器都是稳定的。
而IIR（它是递归滤波器）的极点由分母多项式所确定，极点分布在单位圆内及单位圆上，这意味着不能保证其稳定性，但是从基本信息处可以看到上面所讨论的几种IIR滤波器是稳定的。
4.通常的几种模拟滤波器的特性：
Butterworth滤波器在通带和阻带都是单调的，意味着它们在一个方向上平滑变化；ChebyshevTypeI滤波器在阻带内是单调的，但是在通带内是有波纹的；ChebyshevTypeII滤波器在通带内是单调的，但是在阻带内是有波纹的；Elliptic滤波器在通带和阻带内都有波纹。
5.对于同样的设计要求，切比雪夫滤波器（minorder32）比巴特沃斯滤波器（minorder175）的阶数要低得多。
补充说明：
FIR和IIR比较小结
FIR滤波器
IIR滤波器
H（n）有限长
H（n）无限长
极点固定在原点
极点位于Z平面任意位置
滤波器阶次高
滤波器阶次低
可严格的线性相位
非线性相位
一般采用非递归结构
递归结构
可以用FFT计算
不能用FFT计算
设计借助于计算机
可用模拟滤波器设计
可设计各种幅频特性和相频特性的滤波器
可用于设计规格化的选频滤波器
实验二利用函数法设计FIR及IIR滤波器
一、实验目的
1.熟悉掌握设计滤波器的操作步骤
2.学会用频率响应采样法和窗函数法设计FIR滤波器，并比较其性能
3.学会用单位脉冲响应法和双线性法设计IIR滤波器，并比较其性能
二、实验条件
Matlab2013bpc机
三、实验过程
滤波器指标：
Fs=8000hz；fpass=3000hz；fstop=3050；δs=0.1；δp=0.1
（1）窗函数法设计FIR低通滤波器（主要可以参见实验四）
设计步骤
1.在过渡带宽度的中间，选择通带边缘频率（Hz）：
f1=fpass+（fstop-fpass）/2=3025
2.计算Ω1=2*pi*f1/Fs,并将此值代入理想低通滤波器的脉冲响应h1[n]中：
h1[n]=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中选择满足阻带衰减及其他滤波器要求的窗函数—矩形窗，用表中的公式计算所需窗的非零项目数目
N=0.91*Fs/（fstop-fpass）=145.6
选择奇数项145，这样脉冲响应完全可以对称，避免了滤波器产生相位失真，对于|n|<=（N-1）/2,设计窗函数w[n]。
表1
4.对于|n|<=（N-1）/2,从式h[n]=h1[n]w[n]计算（有限）脉冲响应，对于其他n值h[n]=0,此脉冲响应是非因果的。
5.将脉冲响应右移（N-1）/2，确保第一个非零值在n=0处，使此低通滤波器为因果的。
编写代码如下：
%LPFdesign窗函数法
clc;clearall;closeall;
omga_d=3025*pi/4000;%理想低通滤波器的频率
omga=[0:
pi/100:
pi];
N=145;
w1=window（@rectwin,N）;%矩形窗
subplot（411）;
plot（1:
N,w1）;
axis（[1N01.1]）;%取值范围不同，矩形窗的中心位置不同
title（'窗函数时域形状'）;
legend（'rectwin'）;%标注窗函数名称
M=floor（N/2）;%M=72
n=1:
M;
%hd=sin（n*omga_d）./（n*omga_d）*omga_d/pi;
hd=sin（n*omga_d）./（n*pi）;%理想低通滤波器对应的脉冲响应的右边序列
hdd=[fliplr（hd）,1/omga_d,hd];%绘制原脉冲函数图像，关于（N+1）/2对称
subplot（412）;
plot（1:
N,hdd）;axis（[1N-0.20.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运算得到H
subplot（413）;
plot（omga,H_d1）;gridon;
axis（[0pi01]）;
title（'窗函数法得滤波器形状H'）;
subplot（414）;
plot（abs（fft（h_d1））'）;
title（'fft求得滤波器形状'）;
图形显示：
虽然矩形窗得到的低通滤波器有较好的特性，但是由于其边缘陡峭，旁瓣出现振铃，且变化较为明显。
（2）频率响应采样法设计FIR低通滤波器。
设计原理及滤波器性能分析：
.频率采样法是从频域出发，对给定的理想滤波器的频响
进行N点等间隔采样
然后以此Hd（k）作为实际FIR滤波器的频率特性采样值H（k），即令：
由DFT定义，可以用这N个频域的采样值H（k）来惟一确定FIR的单位脉冲响应h（n），即：
图示如下：
对设计出的滤波器频率响应特性进行分析。
由频域采样定理中的内插公式可以知道，利用这N个频域采样值H（k）同样可以求得FIR滤波器的频率响应
，这个
将逼近理想滤波器的频响。
的内插公式为：
（1）
式中（w）是内插函数：
（2）
把
（2）式带入
（1）式，化简后可得：
（3）
从式（3）可以看到，在各频率采样点上，设计的滤波器，实际的频率响应严格地与理想滤波器的频率响应数值相等，即
但是在采样点之间的频率响应是由各采样点的加权内插函数叠加而形成的，因而有一定的逼近误差。
该误差大小取决于理想频率响应的形状，理想频响特性变化越平缓，内插值越接近理想值，逼近误差越小；反之，如果采样点之间的理想频响特性变化越陡，则内插值与理想值之间的误差越大，因而在理想滤波器不连续点的两边，就会产生尖峰，而在通带和阻带就会产生波纹。
线性相位条件：
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=0~2*pi之间对
进行N点等间隔采样，即
，将w代入式（4）、（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代码中需要用到的函数知识
>>helplinspace
linspaceLinearlyspacedvector.
linspace（X1,X2）generatesarowvectorof100linearly
equallyspacedpointsbetweenX1andX2.
linspace（X1,X2,N）generatesNpointsbetweenX1andX2.
ForN=1,linspacereturnsX2.
编写代码如下：
%LPFdesign频率响应采样法
%待设计滤波器的基本参数
f_p=3000;%所设计滤波器的通带边缘频率
f_s=3050;%阻带边缘频率
Fs=8000;%采样频率%AllinHz
delta_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/2
H=Mag.*exp（-1j*phi）;%|H|*e^j*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）;%相移e^j*omg*M/2,时域采样平移M/2个点
subplot（211）;
plot（h）;%脉冲响应
axis（[1120-0.20.8]）;
title（'脉冲函数hd（n）'）;
subplot（212）;
plot（abs（fft（h）））;%所设计的滤波器的形状
title（'所设计的滤波器形状'）;
图像显示：
增加过渡带抽样点，可加大阻带衰减，但是导致过渡带变宽
增加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用ej*Omga替换得到数字滤波器的频率响应H（Omga）
编写代码如下：
（用要求的参量得不到需要设计的滤波器，这里更换了参数）
%LPFdesign脉冲响应不变法（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=20*log10（abs（H））;
plot（W,H）,gridon;
xlabel（'rad'）;ylabel（'幅度/db'）;
title（'所设计的滤波器形状'）;
axis（[0pi01.1]）;
运行结果
N=3Wn=876.3856
B=1.0e+08*
0006.7311
A=1.0e+08*
0.00000.00000.01546.7311
图形显示：
（4）双线性法设计IIR滤波器
设计步骤
1.确定待求通带边缘频率fpass=3000hz，待求阻带边缘频率fstop=3050hz，待求阻带衰减20db，通带边缘增益-0.91db。
2.用式子Ω=2π*f/Fs把
升级会员 