IIR biquad real-time filter just output noises

Question

I've implemented a IIR filter but it outputs only noise (like FM radio out of tuning) and I can't see what's wrong.

This is the processing function:

void Filter::process(int size, double *in, double *out)
{
    double xb[2] = { 0, 0 };
    double yb[2] = { 0, 0 };

    for (int n = 0; n < size; n++)
    {
        out[n] = (coeffs->a0 * in[n]) + (coeffs->a1 * xb[0])
               + (coeffs->a2 * xb[1]) - (coeffs->b1 * yb[0])
               - (coeffs->b2 * yb[1]);

        xb[1] = xb[0];
        xb[0] = in[n];
        yb[1] = yb[0];
        yb[0] = out[n];
    }
}

in and out are the buffers of sampled data at 48kHz; and size is 1024.

The coefficients are for $ {Fc \over Fs} = {500 \over 48000 } $ and $ Q = 0.7071 $

a0:  0.00107054
a1:  0.00214108
a2:  0.00107054
b1:  -1.99572
b2:  0.953753

Answer 1

I was unable to get a stable filter using your coefficients. You did not state, but am assuming this is a lowpass filter. Using the calculator here, I got:

coeffs = dict(
    a0=0.0010232172047183973,
    a1=0.0020464344094367946,
    a2=0.0010232172047183973,
    b1=-1.9075008174364765,
    b2=0.91159368625535,
)

Filter Code:

I recast to a Transposed Direct Form 2 (Not necessary, but often recommended for floating point) as:

def filter_process(input_data):
    zb = [0., 0.]

    output_data = []
    for data in input_data:

        output_data.append(coeffs['a0'] * data + zb[0])
        zb[0] = zb[1] + coeffs['a1'] * data - coeffs['b1'] * output_data[-1]
        zb[1] = coeffs['a2'] * data - coeffs['b2'] * output_data[-1]

    return output_data

Test Code:

import matplotlib.pyplot as plt
import numpy as np

def plot_f(f):
    Fs = 44000
    sample = 1000
    x = np.arange(sample)
    y = np.sin(2 * np.pi * f * x / Fs)
    out_y = np.array(filter_process(y))
    plt.plot(x, out_y)

for freq in (200, 500, 1000):
    plot_f(freq)

plt.show()

Results:

Update from Comments:

If the input and output buffers are in fact the same buffer, that definitely will be problematic. The transfer function as currently organized, writes a new output to the buffer prior to finishing it use as an input. But that is easily remedied:

void Filter::process(int size, double *in, double *out)
{
    double xb[2] = { 0, 0 };
    double yb[2] = { 0, 0 };

    for (int n = 0; n < size; n++)
    {
        double in_n = in[n]
        out[n] = (coeffs->a0 * in_n) + (coeffs->a1 * xb[0])
               + (coeffs->a2 * xb[1]) - (coeffs->b1 * yb[0])
               - (coeffs->b2 * yb[1]);

        xb[1] = xb[0];
        xb[0] = in_n;
        yb[1] = yb[0];
        yb[0] = out[n];
    }
}

Answer 2

It looks like you got your a and b coefficients swapped. The feedback terms: b0,b1,b2, should be tiny (around 0.00x), and a0,a1,a2 can be close to 1 or 2 in magnitude.

Answer 3

The states should be declared outside the method and be members of your class so they are actually state-full:

class Filter
{
    double xb[2] = { 0, 0 };
    double yb[2] = { 0, 0 };

public:
    void process(int size, double *in, double *out)
    {
        for (int n = 0; n < size; n++)
        {
            out[n] = (coeffs->a0 * in[n]) + (coeffs->a1 * xb[0])
                   + (coeffs->a2 * xb[1]) - (coeffs->b1 * yb[0])
                   - (coeffs->b2 * yb[1]);

            xb[1] = xb[0];
            xb[0] = in[n];
            yb[1] = yb[0];
            yb[0] = out[n];
        }
    }
};