基于FPGA的四阶IIR数字滤波器

【本文由PB创新网为您整理】
摘要：采用FPGA实现四阶IIR数字滤波器，通过两个二阶节级联构成数字椭圆低通滤波器。通带内波纹小于0.1dB，阻带衰减大于32dB。

    关键词：四阶 IIR 椭圆滤波器 补码阵列乘法器

常用的数字滤波器有ＦＩＲ数字滤波器和ＩＩＲ数字滤波器。ＦＩＲ数字滤波器具有精确的线性相位特性，在信号处理方面应用极为广泛，而且可以采用事先设计调试好的ＦＩＲ数字滤波器ＩＰ Ｃｏｒｅ来完成设计，例如Ａｌｔｅｒａ公司提供的针对Ａｌｔｅｒａ系列可编程器件的ＭｅｇａＣｏｒｅ，但是需要向Ａｌｔｅｒａ公司购买或申请试用版。另外，对于相同的设计指标，ＦＩＲ滤波器所要求的阶数比ＩＩＲ滤波器高５～１０倍，成本较高，而且信号的延迟也较大。ＩＩＲ滤波器所要求的阶数不仅比ＦＩＲ滤波器低，而且可以利用模拟滤波器的设计成果，设计工作量相对较小，采用ＦＰＧＡ实现的ＩＩＲ滤波器同样具有多种优越性。

ＩＩＲ滤波器主要有巴特沃斯滤波器、切比雪夫滤波器和椭圆滤波器几种。给出了以上三种滤波器实现同样性能指标所需的阶数及阻带衰减的比较，如表１所示。

表1 三种滤波器的性能比较

原  型	阶  数	阻带衰减/dB
巴特沃斯	6	15
切比雪夫I型	4	25
椭圆函数	3	27

由表１可见，椭圆滤波器给出的设计阶数比前两种低，而且频率特性较好，过渡带较窄，但是椭圆滤波器在通带上的非线性相位响应最明显。本系统选用椭圆函数滤波器进行设计。

１ 原理分析

数字滤波器实际上是一个采用有限精度算法实现的线性非时变离散系统，它的设计步骤为：首先根据实际需要确定其性能指标，再求得系统函数Ｈ（ｚ），最后采用有限精度算法实现。
javascript:window.open(this.src);" style="cursor:pointer;"/>
根据需要，本系统的设计指标为：模拟信号采样频率为２ＭＨｚ，每周期最少采样２０点，即模拟信号的通带边缘频率为ｆｐ＝１００ｋＨｚ，阻带边缘频率ｆｓ＝１ＭＨｚ，通带波动Ｒｐ不大于０．１ｄＢ（通带误差不大于５％），阻带衰减Ａｓ不小于３２ｄＢ。换算为数字域指标为：Ｗｐ＝０．１π，Ｗｓ＝０．２π，Ｒｐ＝０．１ｄＢ，Ａｓ＝３２ｄＢ。系统函数Ｈ（ｚ）的计算采用Ｍａｔｌａｂ软件比较方便，其中有两个现成的函数可以使用：ｅｌｌｉｐｏｒｄ（ｗｐ／ｐｉ，ｗｓ／ｐｉ，Ｒｐ，Ａｓ）函数用来计算数字椭圆滤波器的阶次Ｎ和３ｄＢ截止频率ｗｎ，而ｅｌｌｉｐ（Ｎ，Ｒｐ，Ａｓ，ｗｎ）函数可以求得直接型椭圆ＩＩＲ滤波器的各个系数。通过调用以上两个函数计算得到的系统函数Ｈ（ｚ）为：

javascript:window.open(this.src);" style="cursor:pointer;"/>

这是一个四阶ＩＩＲ系统，Ｍａｔｌａｂ计算出该系统的频率响应如图１所示，可见满足设计要求。

如果采用直接型结构实现，需用的乘法器和延迟单元相对较多，而且分子和分母的系数相差较大，需要较多的二进制位数才能实现相应的精度要求。

如果采用二阶节级联实现，一来各基本节的零点、极点可以很方便地单独进行调整，二来可以降低对二进制数位数的要求。给出了一个直接型结构转为级联型结构的ｄｉｒ２ｃａｓ．ｍ文件，利用该函数求得系统函数的级联表达形式为：

Ｈ（ｚ）＝Ｈ１（ｚ）×Ｈ２（ｚ）＝（0.11-0.1041z -1+0.11z -2）/(1-1.58z -1+0.6469z -2)×(0.2464-0.426z -1+0.2464z -2)/(1-1.7753z -1+0.892z -2)

由上式可以看出，每个二阶节的分子、分母系数差异减少了。值得注意的是，在分配二阶节的增益时，要保证每个节不会发生运算溢出，可以先用Ｍａｔｌａｂ软件分析计算来合理安排各节的增益。经过计算，本文采用第一级分配０.１１，第二级分配０.２４６４,可以保证在要求的输入范围，没有数据溢出发生。
javascript:window.open(this.src);" style="cursor:pointer;"/>
２ 系统实现

将第一个二阶节的系统函数表示为差分方程:

ｙ１(ｎ)＝ａ０ｘ(ｎ)－ａ１ｘ(ｎ－１)＋ａ２ｘ(ｎ)＋ｂ０ｙ(ｎ－１)－ｂ１ｙ(ｎ－２)

＝０．１１ｘ(ｎ)－０.１０４１ｘ(ｎ－１)＋０．１１ｘ(ｎ)＋１.５８ｙ(ｎ－１)－０.６４６９ｙ(ｎ－２)

可以看出，一个二阶节的实现需要五次乘法运算、四次加法运算（采用二进制补码将减法运算变为加法运算）。两个二阶节共需要十次乘法运算。虽然现在已有上千万门的ＦＰＧＡ产品可供选用，但是一般应用时全部采用硬件阵列乘法器毕竟不太合适，而如果采用串行乘法器进行分时复用，其工作速度也不太理想。

本文采用一个折中的方法实现，即乘加单元（ＭＡＣ）的乘法器采用阵列乘法器，而不使用串行乘法器，以提高运算速度。需要注意的是，ＭＡＸ＋ｐｌｕｓⅡ的ＬＰＭ库中乘法运算为无符号数的阵列乘法，所以使用时需要先将两个补码乘数转换为无符号数相乘后，再将乘积转换为补码乘积输出。每个二阶节完成一次运算共需要６个时钟周期，而且需采用各自独立的ＭＡＣ实现两级流水线结构，即每个数据经过两个二阶节输出只需要６个时钟周期。

２．１ 系统原理框图

系统原理框图如图２所示，模拟信号经过ＴＬＣ５５１０转换为００Ｈ～ＦＦＨ的二进制数后，送入四阶ＩＩＲ低通滤波器，处理后输出１０位二进制数送ＡＤ７５２０得到双极性的模拟电压输出。

图3 四阶IIR滤波器的顶层原理图

    ２．２ 顶层ＩＩＲ模块

顶层ＩＩＲ模块如图３所示。主要由一个时序控制模块ＩＩＲＣ、两个ＩＩＲ二阶节模块（ＩＩＲ１和ＩＩＲ２）构成。ＩＩＲ模块设计为１０位二进制补码输入，最高位ａｄ９为补码符号位，次高位ａｄ８用于防止运算时的溢出。可见该ＩＩＲ模块实际可以输入９位二进制补码数，但ＴＬＣ５５１０的输出数据为８位，输入到ＩＩＲ模块时，将ａｄ９和ａｄ８引脚均接地，即输入为正极性电压。

ｃｌｒ输入端为异步清零端，高电平有效。当输入时钟ｃｌｋ为１２ＭＨｚ时，ＩＩＲ模块产生一个频率为２ＭＨｚ的ｃｌｋ＿ａｄ输出时钟提供给ＴＬＣ５５１０。输出数据ｄｏｕｔ为１０位二进制补码。ＩＩＲ１和ＩＩＲ２模块构成级联结构。

２．３ ＩＩＲ１和ＩＩＲ２模块

ＩＩＲ１、ＩＩＲ２模块主要由两个模块构成，一个是数据移位模块，在ＣＬＫ＿Ｒ时钟作用下将差分方程的各ｘ、ｙ值延迟一个时钟；另一个模块是补码乘加单元，用ＶＨＤＬ语言编写，两个乘数先取补后再进行阵列乘法，在ＣＬＫ＿Ｂ时钟控制下完成一次乘加运算，乘积取补后输出，共需要６个时钟。

差分方程的各系数如表２所示，采用１０位定点纯小数补码表示。

表2 二阶差分方程的系数

系  数	a0	a1	a2	b0	b1
IIR1	01CH	3E6H	01CH	194H	35BH
IIR2	03FH	393H	03FH	1C6H	31CH

另外?熏模块中的五个系数定义为常数，以节省硬件资源，并且采用０舍１入法进行数据处理，尽量提高数据运算精度。ＶＨＤＬ程序如下：

ｅｎｔｉｔｙ ｓｍｕｌｔａｄｄ１ ｉｓ

ｐｏｒｔ (ｃｌｋ＿ｒｅｇｂｔ,ｃｌｋ＿ｒｅｇ: ｉｎ ｓｔｄ＿ｌｏｇｉｃ:

ｘ０,ｘ１,ｘ２,ｙ０,ｙ１:ｉｎ ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(９ ｄｏｗｎｔｏ ０);

ｙｏｕｔ: ｏｕｔ ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(９ ｄｏｗｎｔｏ ０));

ｅｎｄ ｓｍｕｌｔａｄｄ１;

ａｒｃｈｉｔｅｃｔｕｒｅ ｂｅｈａｖ ｏｆ ｓｍｕｌｔａｄｄ１ ｉｓ

ｓｉｇｎａｌ ｔａｎ,ｔｂｎ,ｔｐ２ｎ:ｓｔｄ＿ｌｏｇｉｃ;

ｓｉｇｎａｌ ｃｎｔ: ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(２ ｄｏｗｎｔｏ ０);

ｓｉｇｎａｌ ｔａ,ｔｂ,ｔａａ,ｔｂｂ:ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(８ ｄｏｗｎｔｏ ０);

ｓｉｇｎａｌ ｔｍｐａ,ｔｍｐｂ:ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(９ ｄｏｗｎｔｏ ０);

ｓｉｇｎａｌ ｔｐ:ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(１８ ｄｏｗｎｔｏ ０);

ｓｉｇｎａｌ ｔｐｐ:ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ,２２ ｄｏｗｎｔｏ ０);

ｓｉｇｎａｌ ｙｔｍｐ,ｐ:ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(２３ ｄｏｗｎｔｏ ０);

ｃｏｎｓｔａｎｔ ａ０:ｓｔｄ＿ｌｏｇｉｃ＿ｖｅｃｔｏｒ(９ ｄｏｗｎｔｏ ０:＝“０００００１１１００”

（其余常数说明略）

ｂｅｇｉｎ

ｔｐ２ｎ＜＝ｔａｎ ｘｏｒ ｔｂｎ;－－求补后送阵列乘法器

ｔａａ＜＝ｎｏｔ ｔａ ＋‘１’ ｗｈｅｎ (ｔａｎ＝‘１’) ｅｌｓｅ ｔａ;

ｔｂｂ＜＝ｎｏｔ ｔｂ ＋‘１’ ｗｈｅｎ (ｔｂｎ＝‘１’) ｅｌｓｅ ｔｂ;

ｔｐｐ＜＝‘１’＆‘１’＆‘１’＆‘１’＆ ｎｏｔ ｔｐ ＋‘１’ ｗｈｅｎ(ｔｐ２ｎ＝‘１’) ｅｌｓｅ ｔｐ;

ｔｍｐａ＜＝ａ０ ｗｈｅｎ ｃｎｔ＝０ ｅｌｓｅ

ａ１ ｗｈｅｎ ｃｎｔ＝１ ｅｌｓｅ

ａ２ ｗｈｅｎ ｃｎｔ＝２ ｅｌｓｅ

ｂ０ ｗｈｅｎ ｃｎｔ＝３ ｅｌｓｅ

ｂ１ ｗｈｅｎ ｃｎｔ＝４ ｅｌｓｅ (ｏｔｈｅｒｓ＝＞‘０’);

ｔｍｐｂ＜＝ｘ０ ｗｈｅｎ ｃｎｔ＝０ ｅｌｓｅ

ｘ１ ｗｈｅｎ ｃｎｔ＝１ ｅｌｓｅ

ｘ２ ｗｈｅｎ ｃｎｔ＝２ ｅｌｓｅ

ｙ０ ｗｈｅｎ ｃｎｔ＝３ ｅｌｓｅ

ｙ１ ｗｈｅｎ ｃｎｔ＝４ ｅｌｓｅ (ｏｔｈｅｒｓ＝＞‘０’);

ｔａ＜＝ｔｍｐａ(８ ｄｏｗｎｔｏ ０);ｔｂ＜＝ｔｍｐｂ(８ ｄｏｗｎｔｏ ０);

ｔａｎ＜＝ｔｍｐａ(９);ｔｂｎ＜＝ｔｍｐｂ(９);

ｔｐ＜＝ｔａａ*ｔｂｂ;

ｐ＜＝(ｏｔｈｅｒｓ＝＞‘０’) ｗｈｅｎ (ｔｍｐｂ＝“００００００００００”) ｅｌｓｅ

ｔｐ２ｎ ＆ ｔｐｐ;

ｐｒｏｃｅｓｓ (ｃｌｋ＿ｒｅｇ,ｃｌｋ＿ｒｅｇｂｔ)

ｂｅｇｉｎ

ｉｆ ｃｌｋ＿ｒｅｇ＝‘１’ ｔｈｅｎ ｃｎｔ＜＝“０００”;ｙｔｍｐ＜＝(ｏｔｈｅｒｓ＝＞‘０’);

ｅｌｓｉｆ (ｃｌｋ＿ｒｅｇｂｔ’ｅｖｅｎｔ ａｎｄ ｃｌｋ＿ｒｅｇｂｔ＝‘１’) ｔｈｅｎ

ｉｆ ｃｎｔ＜５ ｔｈｅｎ ｃｎｔ＜＝ｃｎｔ＋１;ｙｔｍｐ＜＝ｙｔｍｐ＋ｐ;

ｅｌｓｉｆ (ｃｎｔ＝５) ｔｈｅｎ

ｉｆ ｙｔｍｐ(７)＝‘１’ ｔｈｅｎ

ｙｏｕｔ(８ ｄｏｗｎｔｏ ０)＜＝ｙｔｍｐ(１６ ｄｏｗｎｔｏ ８)＋１;

ｙｏｕｔ(９)＜＝ｙｔｍｐ（２３）；

ｅｌｓｅ ｙｏｕｔ（８ ｄｏｗｎｔｏ ０）＜＝ｙｔｍｐ（１６ ｄｏｗｎｔｏ ８）；

ｙｏｕｔ（９）＜＝ｙｔｍｐ（２３）； ｅｎｄ ｉｆ；

ｅｎｄ ｉｆ；

ｅｎｄ ｉｆ；

ｅｎｄ ｐｒｏｃｅｓｓ；

ｅｎｄ ｂｅｈａｖ；

ＩＩＲ２模块的输出数据采用将补码最高符号位直接取反转换为移码后，就可以送到ＤＡＣ７５２０实现双极性信号输出。
javascript:window.open(this.src);" style="cursor:pointer;"/>
３ 系统性能测试

系统性能的测试采用单极性方波周期信号作为输入信号。信号的频率为１００ｋＨｚ，在采样频率为２ＭＨｚ时，每个周期采样２０个点，换算成数字域频率为０．１π，其二次谐波的数字频率为０．２π。输入到ＴＬＣ５５１０的信号电压幅度为０～２Ｖ，则经过Ａ／Ｄ转换后的输出为００Ｈ～ＦＦＨ。由于低通滤波器的阻带截止频率选在２００ｋＨｚ，衰减３２ｄＢ，由信号理论分析可知，周期方波信号没有二次谐波，所以对三次谐波的衰减经过ＩＩＲ滤波器后输出有直流分量的基波（频率为１００ｋＨｚ）正弦信号。理论计算给出的方波周期信号基波幅度为：

2E/π＝(2×255)/π＝１６２．３４

输入一个周期的数据，Ｍａｔｌａｂ的计算值与ＭＡＸ＋ｐｌｕｓⅡ的仿真值如表３所示。

表3 滤波后输出的数据

输入数据	255	255	255	255	255	255	255	255	255	255
计算值	28.7	-8.2	-29.4	-34.9	-25.2	-1.3	34.8	80.0	130.5	182.0
仿真值	32	1020	999	993	1002	1	36	80	129	179
输入数据	0	0	0	0	0	0	0	0	0	0
计算值	223.4	260.2	281.4	286.9	277.2	253.2	217.1	172.0	121.5	70.1
仿真值	219	255	276	282	273	250	215	171	122	72

由表３可见，仿真输出值为补码，谷点输出值９９３换算成符号数为９９３－１０２４＝－３１。Ｍａｔｌａｂ软件计算的满度输出值为２８６．９，其基波幅度为[２８６.９－(－３４.9)]／２＝１６０.９，与理论值的误差为:

(160.9-162.34)/162.34＝－０.８７％

四阶ＩＩＲ滤波器实现的满度输出值为[２８２－(－３１)]／２＝１５６.５，与理论值的误差为：

（156.5-162.34）/162.34＝－３.６％

这是由于有限精度算法所引起的误差，可以通过增加二进制位数来提高系统的运算精度。图４给出单极性方波信号的前三个周期经过滤波后得到的含直流分量的输出波形，其中实线为Ｍａｔｌａｂ的计算值，“*”为ＭＡＸ＋ｐｌｕｓⅡ的仿真输出。可见，该四阶级联ＩＩＲ滤波器达到了设计要求。

如果改变滤波器的输入时钟频率，则可以改变滤波器的截止频率。另外如果输入无直流分量的周期信号，而且其频率为采样频率的１／２０，则该低通滤波器可以直接得到基波分量输出。其实，要将ＴＬＣ５５１０输出的直流分量滤出很容易，只需利用ＦＰＧＡ做一个减法运算即可。