Implementing a Butterworth Filter Manually in C/C++ via Second Order Sections

Question

Short background: I want to implement a lowpass butterworth filter in C/C++. The end goal is to use this in a low-latency Python program, for which of course scipy.signal exists, but I want to be able to change the cutoff frequency (smoothly) at each time step.

---

Where I am right now: I found a really helpful article on the discretization of a butterworth filter here. It ends with a helpful section on Second Order Sections, where they include the following formula:

And then, to get an n-th order filter, we do $$H_n(s') = \prod_{k=0}^{n/2-1}H_{n,k}(s').$$

So I figured I can just translate the H_s,k into discrete time:

y_k[i] + 2 y_k[i-1] + y_k[i-2] = b0_k x_k[i] + b1_k x_k[i-1] + b2_k x_k[i-2], (Eq. 1)

    where b0_k = gamma^2 - alpha_k gamma + 1,
          b1_k = 2 - 2\gamma^2
          b2_k = gamma^2 + alpha_k gamma + 1,

          gamma, alpha_k as in the above post.

Now comes the tricky part. I actually am a DSP noob, and I have no real clue whether it makes sense to iteratively apply the above. However, I did just that, implementing:

  ... bunch of setup code omitted ...

  for (size_t i = 0; i < num_samples; i++) {
    const float x = input_signal[i];
    const float cutoff_i = cutoffs[i];
    const float gamma = 1/(std::tan(pi*cutoff_i/sampling_rate));
    // We set xk[t] = x[t] for i = {t, t-1, t-2}.
    xk0 = x0;  // x0 is input_signal[i],   xk0 is x_k[i].
    xk1 = x1;  // x1 is input_signal[i-1], xk1 is x_k[i-1].
    xk2 = x2;  // ...
    // Iterate `k` from 0 to `n/2-1`.
    for (size_t k=0; k <= (n/2-1); k++){
      const float alpha_k = 2 * std::cos(2*pi * (2.*k + n + 1)/(4.*n));
      const float b0k = gamma * gamma - alpha_k * gamma + 1.;
      const float b1k = 2. - 2. * gamma*gamma;
      const float b2k = gamma * gamma + alpha_k * gamma + 1.;
      // y0[k] is y_k[i],
      // y1[k] is y_k[i-1],
      // y2[k] is y_k[i-2],
      // Here we calculate y_k[i] according to Eq 1.
      y0[k] = -2.*y1[k] - y2[k] + b0k * xk0 + b1k * xk1 + b2k * xk2;

      // Shift xk.
      xk0 = y0[k]; xk1 = y1[k]; xk2 = y2[k];

      // Shift ys.
      y2[k] = y1[k]; y1[k] = y0[k];
    }

    x2 = x1; x1 = x0;

    output[i] = y0[n/2-1];

---

Sadly, this does not work at all. Outputs blow up to -infinity within a few time steps.

I see a few potential problems: a) The formula for H_2,k above is wrong b) My assumption that I can just go through the sections is wrong c) My code is wrong

Answer 1

c) My code is wrong

That one. You have your difference equations backwards. It should be

$$y[n] = x[n] + 2x[n-1] + x[n-2] - a_1y[n-1] - a_2y[n-2]$$

You have your "a" and "b" coefficient swapped.

You can probably do it this way, but it feels needlessly complicated. IMO it's easier to

1. Start with the poles of an analog prototype filter. These are uniformly distributed on the unit circle on the left side.
2. Pick a sample period that moves the prototype filter to the desired cutoff frequency $ T = 2tan(\pi/2*w) $ where $w$ is the ratio of the desired cutoff to the Nyquist frequency.
3. Transform the analog poles into the z-plane using the bilinear transform: $p_z = \frac{1+p_sT/2}{1-p_sT/2}$
4. Create one section for each conjugate complex pole pair. The forward (or "b") coefficients are simply $b_0 = 1,b_1 = 2, b_2=1$ for all sections. The feedback (or "a" coefficients are $a_0 = 1, a_1 = 2 \cdot \mathfrak{Re}\{p_k\}, a_2 = |p_k|^2$, where $p_k$ is the pole for section k.
5. Normalize the overall gain to 1 at DC

Answer 2

My code is wrong

Even without assuming that the code's behavior is wrong, for long-term maintainability it has its problems.

You'd do much better to structure your code such that you have a data type defined that describes a 2nd-order filter (in C or C++ it would be a struct or class), and a data type that describes a filter's state (in C or C++ it could be an array or struct -- and maybe combined with your filter description if you're only ever going to run on a big machine, but maybe not if you're going to run on a machine that has constant memory).

Make a separate function that executes one 2nd-order filter. Then test the bejabbers out of it, on its own. Test it with all sorts of filters (highpass, lowpass, bandpass, and bandstop) and make sure it works. Only when you're sure it works, then incorporate it into your larger body of code.

So a very C-ish, not very pseudo, pseudo-code of your program would read something like that below.

typedef struct filter_def_t { ... } filter_def_t;
typedef struct filter_state_t { ... } filter_state_t;
double update_filter(
  const filter_def_t * filter,
  filter_state_t * state,
  double input);

  // setup
  filter_def_t filters[N];
  filter_state_t filter_states[N];
  for (int i = 0; i < N; i++) {
    // setup state and filters
  }

  ...

  // Iterate through the filter cascade
  // overall_input is what you want to filter
  double intermediate_result = overall_input;
  for (int i = 0; i < N; i++) {
    intermediate_result = update_filter(
      filters + i,
      filter_states + i,
      intermediate_result);
  }
  // Answer is in 'intermediate_result'
```