基于FPGA的四阶IIR数字滤波器怎么样？

　　常用的数字滤波器有FIR数字滤波器和IIR数字滤波器｡FIR数字滤波器具有精确的线性相位特性，在信号处理方面应用极为广泛，而且可以采用事先设计调试好的FIR数字滤波器IPCore来完成设计，例如Altera公司提供的针对Altera系列可编程器件的MegaCore，但是需要向Altera公司购买或申请试用版｡另外，对于相同的设计指标，FIR滤波器所要求的阶数比IIR滤波器高5~10倍，成本较高，而且信号的延迟也较大｡IIR滤波器所要求的阶数不仅比FIR滤波器低，而且可以利用模拟滤波器的设计成果，设计工作量相对较小，采用FPGA实现的IIR滤波器同样具有多种优越性｡
　　IIR滤波器主要有巴特沃斯滤波器､切比雪夫滤波器和椭圆滤波器几种｡给出了以上三种滤波器实现同样性能指标所需的阶数及阻带衰减的比较，如表1所示｡
　　表1 三种滤波器的性能比较
　　
　　由表1可见，椭圆滤波器给出的设计阶数比前两种低，而且频率特性较好，过渡带较窄，但是椭圆滤波器在通带上的非线性相位响应最明显｡本系统选用椭圆函数滤波器进行设计｡
　　1 原理分析
　　数字滤波器实际上是一个采用有限精度算法实现的线性非时变离散系统，它的设计步骤为：首先根据实际需要确定其性能指标，再求得系统函数H（z），最后采用有限精度算法实现｡
　　
　　根据需要，本系统的设计指标为：模拟信号采样频率为2MHz，每周期最少采样20点，即模拟信号的通带边缘频率为fp=100kHz，阻带边缘频率fs=1MHz，通带波动Rp不大于0.1dB（通带误差不大于5%），阻带衰减As不小于32dB｡换算为数字域指标为：Wp=0.1π，Ws=0.2π，Rp=0.1dB，As=32dB｡系统函数H（z）的计算采用Matlab软件比较方便，其中有两个现成的函数可以使用：ellipord（wp/pi，ws/pi，Rp，As）函数用来计算数字椭圆滤波器的阶次N和3dB截止频率wn，而ellip（N，Rp，As，wn）函数可以求得直接型椭圆IIR滤波器的各个系数｡通过调用以上两个函数计算得到的系统函数H（z）为：
　　
　　这是一个四阶IIR系统，Matlab计算出该系统的频率响应如图1所示，可见满足设计要求｡
　　如果采用直接型结构实现，需用的乘法器和延迟单元相对较多，而且分子和分母的系数相差较大，需要较多的二进制位数才能实现相应的精度要求｡
　　如果采用二阶节级联实现，一来各基本节的零点､极点可以很方便地单独进行调整，二来可以降低对二进制数位数的要求｡给出了一个直接型结构转为级联型结构的dir2cas.m文件，利用该函数求得系统函数的级联表达形式为：
　　H（z）=H1（z）×H2（z）=（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）
　　由上式可以看出，每个二阶节的分子､分母系数差异减少了｡值得注意的是，在分配二阶节的增益时，要保证每个节不会发生运算溢出，可以先用Matlab软件分析计算来合理安排各节的增益｡经过计算，本文采用第一级分配0.11，第二级分配0.2464，可以保证在要求的输入范围，没有数据溢出发生｡
　　
　　2 系统实现
　　将第一个二阶节的系统函数表示为差分方程：
　　y1（n）=a0x（n）-a1x（n-1）+a2x（n）+b0y（n-1）-b1y（n-2）
　　=0.11x（n）-0.1041x（n-1）+0.11x（n）+1.58y（n-1）
　　-0.6469y（n-2）
　　可以看出，一个二阶节的实现需要五次乘法运算､四次加法运算（采用二进制补码将减法运算变为加法运算）｡两个二阶节共需要十次乘法运算｡虽然现在已有上千万门的FPGA产品可供选用，但是一般应用时全部采用硬件阵列乘法器毕竟不太合适，而如果采用串行乘法器进行分时复用，其工作速度也不太理想｡
　　本文采用一个折中的方法实现，即乘加单元（MAC）的乘法器采用阵列乘法器，而不使用串行乘法器，以提高运算速度｡需要注意的是，MAX+plusⅡ的LPM库中乘法运算为无符号数的阵列乘法，所以使用时需要先将两个补码乘数转换为无符号数相乘后，再将乘积转换为补码乘积输出｡每个二阶节完成一次运算共需要6个时钟周期，而且需采用各自独立的MAC实现两级流水线结构，即每个数据经过两个二阶节输出只需要6个时钟周期｡
　　2.1 系统原理框图
　　系统原理框图如图2所示，模拟信号经过TLC5510转换为00H~FFH的二进制数后，送入四阶IIR低通滤波器，处理后输出10位二进制数送AD7520得到双极性的模拟电压输出｡
　　
　　2.2 顶层IIR模块
　　顶层IIR模块如图3所示｡主要由一个时序控制模块IIRC､两个IIR二阶节模块（IIR1和IIR2）构成｡IIR模块设计为10位二进制补码输入，最高位ad9为补码符号位，次高位ad8用于防止运算时的溢出｡可见该IIR模块实际可以输入9位二进制补码数，但TLC5510的输出数据为8位，输入到IIR模块时，将ad9和ad8引脚均接地，即输入为正极性电压｡
　　clr输入端为异步清零端，高电平有效｡当输入时钟clk为12MHz时，IIR模块产生一个频率为2MHz的clk_ad输出时钟提供给TLC5510｡输出数据dout为10位二进制补码｡IIR1和IIR2模块构成级联结构｡
　　2.3 IIR1和IIR2模块
　　IIR1､IIR2模块主要由两个模块构成，一个是数据移位模块，在CLK_R时钟作用下将差分方程的各x､y值延迟一个时钟;另一个模块是补码乘加单元，用VHDL语言编写，两个乘数先取补后再进行阵列乘法，在CLK_B时钟控制下完成一次乘加运算，乘积取补后输出，共需要6个时钟｡ 差分方程的各系数如表2所示，采用10位定点纯小数补码表示｡
　　表2 二阶差分方程的系数
　　
　　另外模块中的五个系数定义为常数，以节省硬件资源，并且采用0舍1入法进行数据处理，尽量提高数据运算精度｡VHDL程序如下：
　　entity smultadd1 is
　　port （clk_regbt，clk_reg： in std_logic：
　　x0，x1，x2，y0，y1:in std_logic_vector（9 downto 0）;
　　yout： out std_logic_vector（9 downto 0））;
　　end smultadd1;
　　architecture behav of smultadd1 is
　　signal tan，tbn，tp2n:std_logic;
　　signal cnt： std_logic_vector（2 downto 0）;
　　signal ta，tb，taa，tbb:std_logic_vector（8 downto 0）;
　　signal tmpa，tmpb:std_logic_vector（9 downto 0）;
　　signal tp:std_logic_vector（18 downto 0）;
　　signal tpp:std_logic_vector，22 downto 0）;
　　signal ytmp，p:std_logic_vector（23 downto 0）;
　　constant a0:std_logic_vector（9 downto0：=“0000011100”
　　（其余常数说明略）
　　begin
　　tp2n《=tan xor tbn;--求补后送阵列乘法器
　　taa《=not ta +‘1‘ when （tan=‘1’） else ta;
　　tbb《=not tb +‘1‘ when （tbn=‘1’） else tb;
　　tpp《=‘1‘&‘1’&‘1‘&‘1’& not tp +‘1‘ when（tp2n=‘1’） elsetp;
　　tmpa《=a0 when cnt=0 else
　　a1 when cnt=1 else
　　a2 when cnt=2 else
　　b0 when cnt=3 else
　　b1 when cnt=4 else （others=》‘0‘）;
　　tmpb《=x0 when cnt=0 else
　　x1 when cnt=1 else
　　x2 when cnt=2 else
　　y0 when cnt=3 else
　　y1 when cnt=4 else （others=》‘0’）;
　　ta《=tmpa（8 downto 0）;tb《=tmpb（8 downto 0）;
　　tan《=tmpa（9）;tbn《=tmpb（9）;
　　tp《=taa*tbb;
　　p《=（others=》‘0‘） when （tmpb=“0000000000”） else
　　tp2n & tpp;
　　process （clk_reg，clk_regbt）
　　begin
　　if clk_reg=‘1’ then cnt《=“000”;ytmp《=（others=》‘0‘）;
　　elsif （clk_regbt’event and clk_regbt=‘1‘） then
　　if cnt《5 then cnt《=cnt+1;ytmp《=ytmp+p;
　　elsif （cnt=5） then
　　if ytmp（7）=‘1’ then
　　yout（8 downto 0）《=ytmp（16 downto 8）+1;
　　yout（9）《=ytmp（23）;
　　else yout（8 downto 0）《=ytmp（16 downto 8）;
　　yout（9）《=ytmp（23）; end if;
　　end if;
　　end if;
　　end process;
　　end behav;
　　IIR2模块的输出数据采用将补码最高符号位直接取反转换为移码后，就可以送到DAC7520实现双极性信号输出｡
　　
　　3 系统性能测试
　　系统性能的测试采用单极性方波周期信号作为输入信号｡信号的频率为100kHz，在采样频率为2MHz时，每个周期采样20个点，换算成数字域频率为0.1π，其二次谐波的数字频率为0.2π｡输入到TLC5510的信号电压幅度为0~2V，则经过A/D转换后的输出为00H~FFH｡由于低通滤波器的阻带截止频率选在200kHz，衰减32dB，由信号理论分析可知，周期方波信号没有二次谐波，所以对三次谐波的衰减经过IIR滤波器后输出有直流分量的基波（频率为100kHz）正弦信号｡理论计算给出的方波周期信号基波幅度为：
　　2E/π= （2×255）/π=162.34
　　输入一个周期的数据，Matlab的计算值与MAX+plusⅡ的仿真值如表3所示｡
　　表3 滤波后输出的数据
　　
　　由表3可见，仿真输出值为补码，谷点输出值993换算成符号数为993-1024=-31｡Matlab软件计算的满度输出值为286.9，其基波幅度为［286.9-（-34.9）］/2=160.9，与理论值的误差为：
　　（160.9-162.34）/162.34=-0.87%
　　四阶IIR滤波器实现的满度输出值为［282-（-31）］/2=156.5，与理论值的误差为：
　　（156.5-162.34）/162.34=-3.6%
　　这是由于有限精度算法所引起的误差，可以通过增加二进制位数来提高系统的运算精度｡图4给出单极性方波信号的前三个周期经过滤波后得到的含直流分量的输出波形，其中实线为Matlab的计算值，“*”为MAX+plusⅡ的仿真输出｡可见，该四阶级联IIR滤波器达到了设计要求｡
　　如果改变滤波器的输入时钟频率，则可以改变滤波器的截止频率｡另外如果输入无直流分量的周期信号，而且其频率为采样频率的1/20，则该低通滤波器可以直接得到基波分量输出｡其实，要将TLC5510输出的直流分量滤出很容易，只需利用FPGA做一个减法运算即可｡