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
| improve this question | | | | |
$\endgroup$
  • 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$ – robert bristow-johnson Mar 13 '17 at 4:14
  • 1
    $\begingroup$ i'm just curious: what does the letter b in your xb[2] and yb[2] states signify? $\endgroup$ – robert bristow-johnson Mar 13 '17 at 4:18
  • $\begingroup$ @robertbristow-johnson In my crazy head the b stands for buffer $\endgroup$ – Victor Aurélio Mar 15 '17 at 5:02
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];
    }
}
| improve this answer | | | | |
$\endgroup$
  • $\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$ – Victor Aurélio Mar 12 '17 at 20:09
  • $\begingroup$ Using coefficients generated from calculator you pointed out and your transfer function and got same result. $\endgroup$ – Victor Aurélio Mar 12 '17 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$ – Stephen Rauch Mar 12 '17 at 20:42
  • $\begingroup$ The input comes from JackAudio to JackAudio. That's, its processing the system audio (music). $\endgroup$ – Victor Aurélio Mar 12 '17 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$ – Stephen Rauch Mar 12 '17 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.

| improve this answer | | | | |
$\endgroup$
  • 2
    $\begingroup$ I cant decide whats a or b some people use a as zeros others uses poles. $\endgroup$ – Victor Aurélio Mar 12 '17 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 '17 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$ – Victor Aurélio Mar 13 '17 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];
        }
    }
};
| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ That didn't help me, same result. $\endgroup$ – Victor Aurélio Mar 14 '17 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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