Practical IIR filter implementation

Question

I'm trying to implement a digital IIR filter on an FPGA and would be happy for some inputs regarding the actual digital implementation. I don't have a lot of experience with FPGAs and digital filters in general and this whole matter seems more difficult than I first expected it to be, especially things like required bitwidths, rounding error/coefficient quantization and saturation. The system uses a Zynq xc7z030 FPGA (speed grade 3) and a clock frequency of 125 MHz (the datarate is 125 MSPS).

I started by looking for a suitable filter design using Matlab's "filter designer". I don't have "hard" specifications (just want to get rid of some noise which occurs around 10 MHz), so I decided to use a 2nd order butterworth filter with a corner frequency of 5 MHz as a first attempt. The resulting coefficients I obtained are:

[b0 b1 b2] = [13.3592e-003 26.7184e-003 13.3592e-003] [a0 a1 a2] = [1.0000e+000 -1.6475e+000 700.8968e-003]

I know that there is different filter structures (Direct form I and II and their transposed versions) to implement this. The problem I face here is timing: I observed that I can only fit two numerical operations (addition or multiplication) into a 8 ns clock period. However, all of those structures usually have more than two additions in series, making it impossible to meet timing (e.g. see here: https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html) and I can't do pipelining in the all pole part since my a1 coefficient is nonzero. I therefore came up with the structure in the attached picture.

I than implemented the filter in VHDL and simulated it. The result is actually quite promising and seems to work fine.

The actual questions I have now are:

• What filter structure would you use for such a filter and why? Is the one I implemented suitable or would you recommend something else? Would it make a difference if I swap the order of the all zero and all pole part?
• My input and output is 16 bits and I used 25 bit to quantize the coefficients because the DSP slices of the Zynq have 25 bit multipliers. For the internal quantities I used 24 bits. The reason for this 24 bits is a little bit arbitrary: I did a bodeplot and a step response of the all pole part and saw that its peak gain is 25.4 dB and the step response has an overshoot up to 19.6. To cover this, my filter internals need 5 bits more than my input -> 21 bits. I added 3 more to be on the save side. How would you properly determine the required bitwidth for the internal quantities? Is there some way to guarantee that no overflow will occur? From the simulations it seems to work, but I'd like to be sure that there is no input signal I'm not aware of which could cause an overflow.
• The coefficients have 23 fractional bits, multiplication is done by multiplying the 25 bits and right shifting the result 23 times. Afterwards, rounding is done by just truncating the last 25 bits. I know that this is not the best way to round the results. How would this be properly done, especially when I have to meet the timing requirements (only one multiplication + one addition allowed in the feedback path)?
• Any other remarks or things I overlook here?

VHDL code:

library ieee;
use ieee.std_logic_1164.all;
use ieee.std_logic_unsigned.all;
use ieee.numeric_std.all;

entity IIR2OrderFilter is
    port (  
        Clk           : in std_logic;         -- Clock
        Rst           : in std_logic;         -- Reset
        InData        : in std_logic_vector(15 downto 0);   -- Data input
        DataOut       : out std_logic_vector(15 downto 0)   -- Data output

    );
end entity IIR2OrderFilter;

architecture behav of IIR2OrderFilter is

Constant C_POLE_COEFF_WIDTH                      : Integer := 25;
Constant C_POLE_SCALE_FACTOR_NUMBIT              : Integer := 23;

Constant C_ZERO_COEFF_WIDTH                      : Integer := 25;
Constant C_ZERO_SCALE_FACTOR_NUMBIT              : Integer := 23;

Constant C_POLE_INTERNAL_BITWIDTH                : Integer := 24;

Constant C_ZERO_INTERNAL_BITWIDTH                : Integer := 24;

Constant CoeffB0                            : Integer := 112065;
Constant CoeffB1                            : Integer := 224130;
Constant CoeffB2                            : Integer := 112065;
-- Constant CoeffA0                         : Integer := ;
Constant CoeffA1                            : Integer := 13819896;
Constant CoeffA2                            : Integer := -5879548;



signal AccuxDN, AccuxDP                     : signed(C_POLE_INTERNAL_BITWIDTH-1 downto 0);
signal AccuPreSum                           : signed(C_POLE_INTERNAL_BITWIDTH-1 downto 0);

signal AccuA1                               : signed(C_POLE_INTERNAL_BITWIDTH-1 downto 0);
signal AccuA2, AccuA2Del                    : signed(C_POLE_INTERNAL_BITWIDTH-1 downto 0);
signal AccuB0                               : signed(C_ZERO_INTERNAL_BITWIDTH-1 downto 0);
signal AccuB1                               : signed(C_ZERO_INTERNAL_BITWIDTH-1 downto 0);
signal AccuB2, AccuB2Del                    : signed(C_ZERO_INTERNAL_BITWIDTH-1 downto 0);

signal AccuB1B2Sum                          : signed(C_ZERO_INTERNAL_BITWIDTH-1 downto 0);
signal AccuB1B2SumDel                       : signed(C_ZERO_INTERNAL_BITWIDTH-1 downto 0);
signal AccuB0B1B2Sum                        : signed(C_ZERO_INTERNAL_BITWIDTH-1 downto 0);


begin


AccuPreSum <= AccuA2Del + signed(InData);
AccuxDN <= AccuPreSum + AccuA1;


IIRRegs: process (Clk)
begin
if rising_edge(Clk) then
    if Rst = '0' then
        AccuxDP <= (others => '0');
        AccuA2Del <= (others => '0');
    else
        AccuxDP <= AccuxDN;
        AccuA2Del <= AccuA2;
    end if;
end if;
end process IIRRegs;

AccuA1 <= resize(shift_right(AccuxDP * to_signed(CoeffA1, C_POLE_COEFF_WIDTH), C_POLE_SCALE_FACTOR_NUMBIT), C_POLE_INTERNAL_BITWIDTH);
AccuA2 <= resize(shift_right(AccuxDP * to_signed(CoeffA2, C_POLE_COEFF_WIDTH), C_POLE_SCALE_FACTOR_NUMBIT), C_POLE_INTERNAL_BITWIDTH);

AccuB0 <= resize(shift_right(AccuxDP * to_signed(CoeffB0, C_ZERO_COEFF_WIDTH), C_ZERO_SCALE_FACTOR_NUMBIT), C_ZERO_INTERNAL_BITWIDTH);
AccuB1 <= resize(shift_right(AccuxDP * to_signed(CoeffB1, C_ZERO_COEFF_WIDTH), C_ZERO_SCALE_FACTOR_NUMBIT), C_ZERO_INTERNAL_BITWIDTH);
AccuB2 <= resize(shift_right(AccuxDP * to_signed(CoeffB2, C_ZERO_COEFF_WIDTH), C_ZERO_SCALE_FACTOR_NUMBIT), C_ZERO_INTERNAL_BITWIDTH);



AccuB1B2Sum <= AccuB2Del + AccuB1;
AccuB0B1B2Sum <= AccuB1B2SumDel + AccuB0;


FIRRegs: process(Clk)
begin
if rising_edge(Clk) then
    if Rst = '0' then
        AccuB2Del <= (others => '0');
        AccuB1B2SumDel <= (others => '0');
        DataOut <= (others => '0');
    else
        AccuB2Del <= AccuB2;
        AccuB1B2SumDel <= AccuB1B2Sum;
        DataOut <= std_logic_vector(resize(AccuB0B1B2Sum, 16));
    end if;
end if;
end process FIRRegs;



end architecture behav;

Answer 1

Is the one I implemented suitable or would you recommend something else?

I think you did pretty well there. I'm not sure I'd end up with something different (although I would do a literature search on implementing IIR filters in FPGAs).

Your filter looks pretty close to a direct form 2 filter. You can look up references to that for its strengths and weaknesses. If you're not already immersed in DSP, I'd suggest "Understanding Digital Signal Processing" by Rick Lyons, but only because I know the book.

Swapping the zero part with the pole part would make it a transposed direct-form 2 filter, which will be discussed in those same references.

How would you properly determine the required bitwidth for the internal quantities?

It looks like you found, by experimentation, the correct data path widths. The rule of thumb for setting the data path widths of IIR filters is, for each filter section, compute the product of the distances of each pole from the unit circle, then figure out how many extra bits you need, over the incoming data width, to "hold" the extra precision you're gaining. In your case you have poles at approximately $z \simeq 0.824 \pm 0.149j$. These have absolute values of $|z| \simeq 0.837$; $1 - 0.837 \simeq 0.163 \simeq 2^{-2.6}$. Since you have two of them, you need five extra bits -- so the minimum amount of pad you'd want over your 16 bit data would give you 21-bit data paths; 25 is plenty.

Is there some way to guarantee that no overflow will occur?

Yes. Dust off your block diagram computing skills, and calculate the transfer function from the input signal to the output of each summing block (and multiplier, if it's greater than 1). The peak value of the resulting Bode plot is the worst-case gain from input to that node; scale your variables so that with the maximum peak-peak input at exactly that frequency, that node does not overflow.

Afterwards, rounding is done by just truncating the last 25 bits. I know that this is not the best way to round the results. How would this be properly done,

If your output has excess bits, truncation may well be good enough. Rounding is easy to implement, and should preserve the correct average value at least.

For a given system, I'd figure this out by looking at the anticipated noise in the input signal, and comparing it to the quantization noise caused by truncation or rounding. If the signal's inherent noise is four or more times greater (i.e., around 10dB greater), then it'll swamp out the effect of quantization and you're done.

If you find that you're losing valuable data at this point, you may need to think about wider data paths leading out of the filter.

Answer 2

i dunno VHDL but i suspect that this filter is being implemented in fixed-point arithmetic, is that suspicion correct? if so, can you maintain fixed-point words of different word widths? when you multiply an $N_c$-bit coefficient against an $N_s$-bit signal or state, the resulting product is an $N_c \times N_s - 1$ bit word.

if so, are you aware of simple techniques to mitigate quantization error? first, you should do the Direct Form 1 (which swaps your feeback section (with the $a_n$ coefficients) and the feedforward section (with the $b_n$ coefficients). you will have one big summing node where your sum is of five double wide words. the result of that sum is the only place where you must reduce the word width from $N_c \times N_s + 2$ bits (because of addition, the word width might need to be a couple bits bigger) to an $N_s$-bit word.

the other important technique to do really good fixed-point IIR filtering is the use of noise-shaping in the quantization. the simplest form is 1st-order noise shaping with a zero at DC ($z=1$) which makes your filter perfect at DC (no quantization error at all at DC) but more noisy at Nyquist (where $z=-1$).

Answer 3

There is a technique called the scattered look-ahead transformation to transform a second-order IIR filter in an order-4 IIR filter in order to allow pipelining in the feedback path.

The gist of it is to change $H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 + a_1z^{-1} + a_2z^{-2}}$

To this

$H_1(z) = \frac{d_0 + d_1z^{-1} + d_2z^{-2}+d_3z^{-3}+d_4z^{-4}}{1 + c_2z^{-2}+c_4z^{-4}} = H(z)$

since the odd feedback coefficients are 0, it is possible to use pipelining in the feedback loop thus increasing the maximum frequency.

To go from $H(z)$ to $H_1(z)$ ou need to factor $$H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 + a_1z^{-1} + a_2z^{-2}} = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 -2rcos(\theta)z^{-1} + r^2z^{-2}} $$

$$H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{1 -2rcos(\theta)z^{-1} + r^2z^{-2}} * \frac{re^{±(j(\theta+\pi/2))}}{re^{±(j(\theta+\pi/2))}} = \frac{d_0 + d_1z^{-1} + d_2z^{-2}+d_3z^{-3}+d_4z^{-4}}{1 + c_2z^{-2}+c_4z^{-4}}$$

If your poles are not complex, the technique I showed will not work. However you can still apply the look-ahead. You simply need to transform the order-2 filter in 2 cascaded order1- filters and then apply the lookahead transform.

You should consult this site for more information http://people.ece.umn.edu/users/parhi/SLIDES/chap10.pdf

In the second image, the original IIR filter has been transformed into an order-4 IIR filter using the scattered look-ahead technique. Delays (or flip-flops) have been added in the FIR part of the filter to help reach timing closure. In the feedback path, the multipliers are pipelined. It is possible to do so because the $\ c_1 $ and $\ c_3 $ are equal 0.

The transfer function is the same except for a pipeline delay of 6 cycles. Something that's easily manageable to if you use a data valid flag.

Edit : It would be possible to decrease the pipeline delay from 6 cycles to something like 2-3 cycles depending on the FPGA architecture.