2
$\begingroup$

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
$\endgroup$
3
  • 1
    $\begingroup$ you b1 and b2 coefs don't look like they result in stable poles. |b1| needs to be smaller than 1+b2 . |b1+b2| should be less than 1. $\endgroup$ Mar 13, 2017 at 4:14
  • 1
    $\begingroup$ i'm just curious: what does the letter b in your xb[2] and yb[2] states signify? $\endgroup$ Mar 13, 2017 at 4:18
  • $\begingroup$ @robertbristow-johnson In my crazy head the b stands for buffer $\endgroup$ Mar 15, 2017 at 5:02

3 Answers 3

2
$\begingroup$

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:

enter image description here

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];
    }
}
$\endgroup$
13
  • $\begingroup$ Yep, low pass filter, I have used formulas on the EQ Cookbook[0] to generate the coefficients [0]: musicdsp.org/files/Audio-EQ-Cookbook.txt $\endgroup$ Mar 12, 2017 at 20:09
  • $\begingroup$ Using coefficients generated from calculator you pointed out and your transfer function and got same result. $\endgroup$ Mar 12, 2017 at 20:36
  • $\begingroup$ I do not have a C compiler handy, but I hand translated your C to python, and after I changed coefficients, your transfer function worked for me. That leaves as a key difference the input/output data. You did not specify what you were using. $\endgroup$ Mar 12, 2017 at 20:42
  • $\begingroup$ The input comes from JackAudio to JackAudio. That's, its processing the system audio (music). $\endgroup$ Mar 12, 2017 at 20:46
  • $\begingroup$ Would suggest some test data until you get this running, or an all-pass filter to verify your data pipeline. $\endgroup$ Mar 12, 2017 at 20:47
2
$\begingroup$

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.

$\endgroup$
3
  • 2
    $\begingroup$ I cant decide whats a or b some people use a as zeros others uses poles. $\endgroup$ Mar 12, 2017 at 20:07
  • $\begingroup$ @VictorAurélio The feedback terms or poles usually need to be small to keep thing from eventually blowing up (becoming NaNs or Inf). $\endgroup$
    – hotpaw2
    Mar 12, 2017 at 22:52
  • $\begingroup$ What I mean is the variable name, I find some to use b as poles, other people use b as zeros... $\endgroup$ Mar 13, 2017 at 16:43
1
$\begingroup$

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];
        }
    }
};
$\endgroup$
1
  • $\begingroup$ That didn't help me, same result. $\endgroup$ Mar 14, 2017 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.