1、FIR 数字滤波器分布式算法的原理及 FPGA 实现摘要：在利用 FPGA 实现数字信号处理方面，分布式算法发挥着关键作用，与传统的乘积-积结构相比，具有并行处理的高效性特点。详细研究了基于 FPGA、采用分布式算法实现FIR 数字滤波器的原理和方法，并通过 Xilinx ISE 在 Modelsim 下进行了仿真。 关键词：分布式算法 DALUT FPGA FIR数字滤波器正在迅速地代替传统的由 R、L、C 元件和运算放大器组成的模块滤波器并且日益成为 DSP 的一种主要处理环节。FPGA 也在逐渐取代 ASIC 和 PDSP，用作前端数字信号处理的运算（如：FIR 滤波、CORDIC 算法
2、或 FFT）。乘累加运算是实现大多数 DSP 算法的重要途径，而分布式算法则能够大大提高乘累加运算的效能。1 传统的乘累加结构（MAC Multiplying Accumulator）FIR 数字滤波器基本理论FIR 滤波器被称为有限长脉冲响应滤波器，与 IIR 数字滤波器相对应，它的单位脉冲响应h(n)只有有限个数据点。输入信号经过线性时不变系系统输出的过程是一个输入信号与单位脉冲响应进行线性卷积的过程，即：式中，x(n)是输入信号，y(n)是卷积输出，h(n)是系统的单位脉冲响应。可以看出，每次采样 y(n)需要进行 L 次乘法和 L-1 次加法操作实现乘累加之和，其中 L 是滤波器单位脉
3、冲响应 h(n)的长度。可以发现，当 L 很大时，每计算一个点，则需要很长的延迟时间。2 乘累加运算的位宽分配DSP 算法最主要的就是进行乘累加运算。假设采样信号的位宽用 N 来表示，则 N 位与 N 位的乘累结果需要 2N 位的寄存器来保存；如果两个操作数都是有符号数，则乘积只有 2N-1个有效位，因为产生了两个符号位。为了使累加器的结果不产生溢出(计算结果超出了范围，产生了一个进位)，需要对累加器进行冗余设计，也就是说要在累加器 2N 的位宽上多设计出 K 位，累加器的长度 M 计算方式如下（L 为滤波器的长度）：对于无符号数：M=2N+K=2N+log2 L对于有符号数：M=2NK=2N
4、+log2 L-13 乘累加运算的分布式算法原理分析得益于 Xilinx FPGA 查找表结构的潜能，分布式算法在滤波器设计方面显示出了很高的效率，自 20 世纪 90 年代初以来越来越受到人们的重要。分布式算法是基于查找表的一种计算方法，在利用 FPGA 实现数字信号处理方面发挥着重要的作用，可以大大提高信号的处理效率。它主要应用于数字滤波、频率转换等数字信号处理的乘累加运算。分布式算法推导如下：设 Ak 是已知常数（如滤波器系数 Taps、FFT 中的正弦/余弦基本函数等），Xk(n)是变量，可以看作是 n 时刻的第 k 个采样输入数据，y(n)代表 n 时刻的系统响应。那么它们的内积为：
5、其中，Xk(n)变量可以写成下面的格式： 式中，B 为数据格式的字长, Xkb是变量的二进制位，只有“0”和“1”两种状态。将（2）式代入（1）式得：（这种方式把 MAC 运算转化为根据输入信号 x 寻址 ROM 的加法运算，需要 2N 的 ROM 空间,N 为系数的个数，和 2n 相乘只需左右移位，系数 Ak 如果是负数，就使用补码表示，不对，减法怎么做？）4 FPGA 实现过程中查找表的构造方法根据以上论述，括号中的每一乘积项代表着输入变量的某一位与常量的二进制“与”操作，加号代表着算术和操作，指数因子对括号中的值加权。如果事先构造一个查找表，该表存储着括号中所有可能的组合值，就可以通过所
6、有输入变量相对应位的组合向量（XNb,X(N-1)b,.x1b）对该表进行寻址，该查找表称为 DALUT。DALUT 的构造规则如表 1 所示。5 采用分布式算法实现 FIR 数字滤波器为了说明问题，以一个三个系数的 FIR 数字滤波器为例设计分布式算法，字宽也设置为三位。设 FIR 数字滤波器系数为：h(0)=5,h(1)=2,h(2)=3。在进行 FPGA 设计时，该表以组件 Component 形式构建，设置为 ROM 结构，提供输入寻址端口 table_in2.0,输出端口 table_out3.0。FPGA 算法的结构图如图 2 所示。算法实现中的几个关键问题为：（1）采用状态机实现
7、分布式算法的状态转移状态机的实现如图 3 所示，设置三个状态 s0、s1、s2 。状态 s0 完成数据的装入，数据寄存器需要成对出现，一个完成数据的延迟，另一个完成数据的移位，并将状态转移到 s1；状态 s1 完成查找表功能、数据移位和分布式算法的乘累加运算，数据移位一个数据宽带后将状态转移到 s2;状态 s2 完成数据的输出，并将状态转移到 s0。利用状态机可以条理清楚地简化计算过程，在算法实现时发挥着关键的作用。（2）系统时钟与数据输入时钟的关系根据上述的状态转移关系，可以得出：每输入一个数据，在下一次数据输入之前，需要在状态 s1 停留一个数据宽带（三位）的时钟时间，在 s2 停留一个时
8、钟的数据输出时间。也就是说，系统时钟频率应是数据输入频率的 5 倍，即 fclkock=5fxin。（3）分布式算法中的乘累加式公推导及核心代表实现设 B 是数据的字宽，Pn 是分布式算法第 n 位的结果，则有：有了该关系式，就可以通过 for.loop 循环，使用一条语句完成分布式乘累加算法。具体如下：for n in 0 to B-1 loopP:=p/2+tableout(n)*2B-1;End loop;6 算法仿真验证与结论本文实现的 FIR 滤波器在 Xilinx 的集成开发环境 ISE 下利用 ModelSim 进行了仿真。当输入数据为 7，3，1.时，仿真输出依次为 35,29
9、,32,16.,与乘累加方式 FIR 滤波算法得出的结果完全一致。假设查找表和 PDSP 的通用乘法器延时时间相同，分布式算法的等待时间是 Br,通用乘法器的等待时间是 N1。可见，对于位宽较小的数据来说，分布式算法的执行速度远高于乘累加运算。可见，利用 FPGA 实现分布式计算大大提高了计算的速度，在高速信号处理中发挥着重要作用。*FIR 中，h(k)系数为 N 个， ，求 y(n) 需要前面的 N 个 x 数据。10)()(Nkknxhny )1(*)1(.2*)(0)( nhxhxny（for example y(0)=h(0)*x(-1)+ h(1)*x(-2) + h(2)*x(-3
10、) N=3）x 有 M 位010.1 .2)(x MMM（for example x(-1) = 1100 = 1*23+1*22+0*2+0） )1(*)1(.)2(*)(*)()( Nnhnxhnxhxny .20 011hMM )2(.)()(*)1 00xxx.+ )2*1(.*)1(2*)1()( 001 NnxNnNnh MMM= xhxhx )()(.)()(*)0 1111 2* MMMM nnnh.+0000 )()(.)()(*) Nxhxhx只要 h 的值确定了，把输入数据 x 的各个位提取出来，做一个 LUT，储存 h(0)+h(N-1)的不同值，需要 2N 个 reg
11、ister。DA 大大简化了运算。ModelSim shortcut keyF2 在各个窗口之间切换F3 把选中的窗口放大F9 运行 run在 wave 窗口察看波形时候的快捷键F4 ：是 zoom full 的快捷键F5 ：是放大快捷键F6 ：是缩小快捷键F7 ：是 zoom last 的快捷键还可以按住 ctrl 用鼠标的从左到右直拉，从右到左直拉，斜拉来实现缩放F10 ：run -continueF11 ：run -stepF12 ：run -over在 wave 里面双击想要查看的波形可以调出其 dataflow 对话框的 show drivers波形比较的时候移动看波形差别的快捷键
12、tab shift +tab在 dataflow 窗口中移动整个设计的方法：ctrl 和鼠标的中间键ctrl+f 查找的快捷键F8 在命令窗口匹配前面输入过的命令F3 查找下一个FFT输入数据 x 24bit，全是无符号数。转换为 float： 32 3124 2300 127+24 xIEEE 的标准 32 bits floating 规定尾数默认小数点前面为 0，我们这里默认小数点前面为0 即可。而且考虑一下指数是否需要 8 位来表示？因为乘法后尾数的精度可能会不够？module FIR_Integer_Taps(x0,reset,clk,h_0,h_1,h_2,y);RAM1RAM2Re
13、ad ADFFTMultiplexer Multiplexerparameter BitSize=7;parameter TapNumber=3; /3 coefficienceinput unsigned BitSize:0 x0; / input x is 8 bits, unsigned, if its signed,must take/ the negative into account.reg unsigned BitSize:0 x_1,x_2,x_3;input signed BitSize:0 h_0,h_1,h_2; /3 taps,not symmetricalinput
14、reset,clk;/output ready; /ready=1 calculation ends, input can change;/ =0 calcalation is processing output signed 2*BitSize+4:0 y;reg signed BitSize:0 h0,h1,h2;reg signed 2*BitSize+4:0 sum,y;reg 1:0 ss;reg 2:0 i;reg 1:0 j;reg TapNumber-1:0 temp;reg signed 7:0 Sum_h 7:0; / array for taps summary LUT/
15、reg 2:0 M; /input data is 7 bit, so M is 3 bit (23=8)always (posedge clk) /Main FSMbeginif (reset=1)begin x_1=0;x_2=0;x_3=0;y=0;i=0;j=0;ss=2b00; h0=h_0;h1=h_1;h2=h_2;sum=0;endelse begincase(ss) /State00: LUT Sum_h calculation 2b00: beginSum_h0=0; /i i2 i1 i0 Sum_h1=h2; / h0 h1 h2Sum_h2=h1;Sum_h3=h1+
16、h2;Sum_h4=h0;Sum_h5=h0+h2; Sum_h6=h0+h1;Sum_h7=h0+h1+h2;ss=01;end2b01: /State01: Read 3 value for x()begin x_3=x0;x_2=x_3;x_1=x_2;if(j2b10) j=j+1;else if(j=2b10) ss=2b11; /temp=x_1i,x_2i,x_3i;end2b11: /State11: sum calculationbegin /8 periods needed to get the result. 输入数据为 L 位，则此处需要 L 个周期。/假设此 modu
17、le 频率 f，则前方的采样数据频率必须低于 f/L。否则应该考虑把此步使用并/行方式实现，以面积资源换速度。temp=x_1i,x_2i,x_3i; if (iBitSize)begincase(temp)3b000: sum=sum; 3b001: sum=sum+(Sum_h1i); /left shift i bit3b010: sum=sum+(Sum_h2i);3b011: sum=sum+(Sum_h3i);3b100: sum=sum+(Sum_h4i);3b101: sum=sum+(Sum_h5i); 3b110: sum=sum+(Sum_h6i);3b111: sum=
18、sum+(Sum_h7i); default:sum=0;endcasei=i+1; endelse beginss=2b10; i=0;end end 2b10: /State10:output the result,back to state01beginy=sum;ss=2b01;sum=0;enddefault: ss=2b00;endcaseendendendmoduleThroughput (Sample Rate) Comparison of Single-MAC-Based FIR and DA FIR as a Function of Filter Length. B is
19、the DA FIR Input Sample Precision. The Clock Rate is 120 MHz.A comparison between a DA FIR architecture and a conventional scheduled MAC-based approach. The clock rate is assumed to be 120 MHz for both filter architectures. Several values of input sample precision for the DA FIR are presented. The d
20、ependency of the DA filter throughput on the sample precision is apparent from the plots. For 8-bit precision input samples, the DA FIR maintains a higher throughput for filter lengths greater than 8 taps. When the sample precision is increased to 16 bits, the crossover point is 16 taps.可以看到 DA 算法占据了明显的优势。DA 支持的 sample rate 仅仅与 sample data 的字长(B)有关，而与 filter length 无关。但是如果 sample data 的一个数据长度 B 太大，其需要的ROM 也是指数增长 2B。So there is a tradeoff between DA and MAC algorithm.
