Help with my first IIR biquad (DF1)

[Dmeads]

Hello, I am reading this book: https://link.springer.com/book/10.1007/978-3-642-45309-0

A few weeks ago, I made a first order IIR (lossy integrator) in verilog. I used the book example, then made my own using fixed point. It worked perfect and simulation results matched the books.

I am moving on to biquads, which the book does not have example code for, and I am really struggling with why my filter is unstable. You can see the output is all over the place in the picture below. I am thinking maybe this is due to roundoff error? or something in the HDL that isnt correct? I am trying to make a Direct-form I biquad, because I have read those are the least sensitive to rounding error. I have included a picture of the structure from the book I am trying to code.



I generated the coefficients in Matlab (script included in the zip), then multplied by 2^14 to get the integer coefficients. My iir filter uses signed 16-bit by signed 16-bit multipliers.

Is there anything glaring in my code?

Code: [Select]

`timescale 1ns / 1ps
/* My first try at implementing an IIR filter on FPGA
*  lowpass elliptical filter
*  Cuttoff frequency of 60 kHz, with 40 dB stopband attenuation and 0.5 dB ripple in the passband
*  Direct-form I structure BiQuad
* 
*  Source 1: https://ccrma.stanford.edu/~jos/filters/BiQuad_Section.html#:~:text=The%20term%20%60%60biquad%27%27%20is%20short%20for%20%60%60bi-quadratic%27%27%2C%20and,be%20called%20the%20overall%20gain%20of%20the%20biquad.
*
*  v1.0 10/23/2024 Dominic Meads
*/

module iir_DF1_Biquad (
  input clk,
  input signed [15:0] din,
  output signed [15:0] dout
);

  // filter coefficients (multiplied floating point coefficients by 2^14)
  // sos = {1.0000   -1.8057    1.0000    1.0000   -1.9459    0.9480}
  //         b0         b1        b2        a0        a1        a2
  //   g = 0.0102 
  reg signed [15:0] a1_fixed = -31881;
  reg signed [15:0] a2_fixed = 15531;
  reg signed [15:0] b0_fixed = 167;    // g * b0 * 2^14  (multiply denom coeffs by gain for DF1 [source 1])
  reg signed [15:0] b1_fixed = -302;   // g * b1 * 2^14 
  reg signed [15:0] b2_fixed = 167;    // g * b2 * 2^14

  // input register
  reg signed [15:0] r_din = 0;

  // output register
  reg signed [15:0] r_dout = 0;

  // delay registers (after multiplication)
  reg signed [15:0] z2 = 0;
  reg signed [15:0] z1 = 0;

  // multiplication wires
  wire signed [31:0] w_product_a1;
  wire signed [31:0] w_product_a2;
  wire signed [31:0] w_product_b0;
  wire signed [31:0] w_product_b1;
  wire signed [31:0] w_product_b2;

  // truncated (divided) multiplication wires
  wire signed [15:0] w_trunc_product_a1;
  wire signed [15:0] w_trunc_product_a2;
  wire signed [15:0] w_trunc_product_b0;
  wire signed [15:0] w_trunc_product_b1;
  wire signed [15:0] w_trunc_product_b2;

  always @ (posedge clk)
    begin
      r_din <= din;
      z1 <= w_trunc_product_a1 + w_trunc_product_b1 + z2; 
      z2 <= w_trunc_product_a2 + w_trunc_product_b2;
      r_dout <= w_trunc_product_b0 + z1;
    end

  assign dout = r_dout;

  // multiply
  assign w_product_a1 = r_dout * a1_fixed;
  assign w_product_a2 = r_dout * a2_fixed;
  assign w_product_b0 = r_din * b0_fixed;
  assign w_product_b1 = r_din * b1_fixed;
  assign w_product_b2 = r_din * b2_fixed;

  // divide multiplications by 2^14
  assign w_trunc_product_a1 = w_product_a1 >>> 14;
  assign w_trunc_product_a2 = w_product_a2 >>> 14;
  assign w_trunc_product_b0 = w_product_b0 >>> 14;
  assign w_trunc_product_b1 = w_product_b1 >>> 14;
  assign w_trunc_product_b2 = w_product_b2 >>> 14;

endmodule

[Dmeads]

Okay, after some good research and diving into simulation/matlab, it seems like I have quite a few bad things going on.

The number one issue is overflow happening during the summing, as well as some code bugs like sizing of a few registers.

The coefficients of the filter, even when rounded, are stable. It is the overflow that is making things act weird.

Additionally, the filter structure picture in the book is not Direct-Form 1, but rather Direct-Form 1 Transposed. In this forum, it is stated that Direct-form 1 Transposed filters "should never be used in fixed-point implementation because the poles come first resulting in overflow in these two cases."

I am sure you can use them, but you need some logic to handle overflows.

There is also this online book https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html That talks about how you can ignore overflow as long as the output signal is in range, but it does not state this in the section on DF1-T, so I am going to assume you cannot ignore overflow for DF1-T

I will work on a DF1 code and post it if it works.

[Someone]

it is stated that Direct-form 1 Transposed filters "should never be used in fixed-point implementation because the poles come first resulting in overflow in these two cases."

That is true only if you use a register size smaller than the largest value that could exist. Various solutions exist to resolve that:
Simple filters can have their internal signal magnitudes calculated or estimated accurately enough to size the registers (generally individually) for no possibility of overflow.
For larger order filters it becomes impractical to put a bound on the internal signal magnitudes, saturation can be added to approximate sizing with little effect on performance.

[Dmeads]

For larger order filters it becomes impractical to put a bound on the internal signal magnitudes, saturation can be added to approximate sizing with little effect on performance.

Yes, I read about saturation being used in a white paper but wasn't sure exactly what it was. Will have to do some more research on it. Thanks for clarifying

[Dmeads]

I SOLVED IT!

After rewriting in DF1 form, i was still getting weird results. It was because I was registering the output in a way that added an extra delay to the recursive side of the filter. I fixed it and everything is good matches MATLAB comparison.

Here is the code:

Code: [Select]

`timescale 1ns / 1ps
/* My first try at implementing an IIR filter on FPGA
*  lowpass elliptical filter
*  Cuttoff frequency of 60 kHz, with 40 dB stopband attenuation and 0.5 dB ripple in the passband
*  Direct-form I structure BiQuad
* 
*  Source 1: https://ccrma.stanford.edu/~jos/filters/BiQuad_Section.html#:~:text=The%20term%20%60%60biquad%27%27%20is%20short%20for%20%60%60bi-quadratic%27%27%2C%20and,be%20called%20the%20overall%20gain%20of%20the%20biquad.
*
*  v1.0 10/23/2024 Dominic Meads
*
*/

module iir_DF1_Biquad (
  input clk,
  input signed [15:0] din,
  output signed [15:0] dout
);

  // filter coefficients (multiplied floating point coefficients by 2^14)
  // sos = {1.0000   -1.8057    1.0000    1.0000   -1.9459    0.9480}
  //         b0         b1        b2        a0        a1        a2
  //   g = 0.0102 
  reg signed [15:0] a1_fixed = -31881;
  reg signed [15:0] a2_fixed = 15531;
  reg signed [15:0] b0_fixed = 167;    // g * b0 * 2^14  (multiply denom coeffs by gain for DF1 [source 1])
  reg signed [15:0] b1_fixed = -302;   // g * b1 * 2^14 
  reg signed [15:0] b2_fixed = 167;    // g * b2 * 2^14

  // input register
  reg signed [15:0] r_x = 0;

  // output register
  reg signed [15:0] r_y = 0;

  // delay registers
  reg signed [15:0] r_x_z1 = 0;
  reg signed [15:0] r_x_z2 = 0;
  reg signed [15:0] r_y_z1 = 0;
  reg signed [15:0] r_y_z2 = 0;

  // multiplication wires
  wire signed [31:0] w_product_a1;
  wire signed [31:0] w_product_a2;
  wire signed [31:0] w_product_b0;
  wire signed [31:0] w_product_b1;
  wire signed [31:0] w_product_b2;

  wire signed [31:0] w_sum; 

  always @ (posedge clk)
    begin
      r_x <= din;
      r_x_z1 <= r_x;
      r_x_z2 <= r_x_z1;
      r_y_z1 <= w_sum >>> 14;  // divide by the same 2^14 value the coefficients were multiplied by
      r_y_z2 <= r_y_z1;
    end

  // multiply
  assign w_product_a1 = r_y_z1 * -a1_fixed;
  assign w_product_a2 = r_y_z2 * -a2_fixed;
  assign w_product_b0 = r_x * b0_fixed;
  assign w_product_b1 = r_x_z1 * b1_fixed;
  assign w_product_b2 = r_x_z2 * b2_fixed;

  assign w_sum = w_product_b0 + w_product_b1 + w_product_b2 + w_product_a1 + w_product_a2;

  assign dout = r_y_z1;

endmodule

[SiliconWizard]

Nice.

[BrianHG]

Q: Is that DC offset on your output intentional?

[Dmeads]

Yes, thank you for checking. Im simulating samples coming from my unipolar ADC on the input. Its a highpass filter, so the DC offset should remain.

The cursor in the waveform screenshot (offscreen) is also in a poor place, because the output is still in the oscillation, while the input is at a DC value, therefore the values look odd on the screenshot. Everything does match up

[glenenglish]

recursive fixed point filters on FPGAs are things are the devil

 avoid  any wraps of course, overflows... Because the whole thing goes to pot even with a clip / saturating arithmetic, do the worst case scenarios for bit width requirement.

and after you conquer all that, you will encounter limit cycles,  LOL which will degrade your SNR, and you'll add another couple of bits precision to the filter to mitigate that (the easy way, without making it a research project). so whatever dynamic range you want, add a couple of bits to it for implementation.

*** there are several tools and filter topologies and strategies  to minimize word sizes / coefficient quantization issues  etc  can help quite a lot. There are books about fixed point recursive filter implementation that explore methods... worthwhile tracking them down.
-glen