0
$\begingroup$

I would like to know the theory behind matlab butter function. Is there any book which gives clear idea on how it works. In future I would like to calculate the coefficients manually.

$\endgroup$
  • $\begingroup$ it would be painful, but i might be able to write a little primer on the theory of the Butterworth filter. maybe limit it to low-pass. Butterworth filters are "maximally flat" in both passband and stopband. this means that in the magnitude frequency response as many derivatives are zero at a frequency of zero (DC). this stuff is in textbooks and maybe some online resource. i wonder if Julius Smith has a good detailed treatment. $\endgroup$ – robert bristow-johnson Apr 3 '15 at 0:37
2
$\begingroup$

The MATLAB function that you refer to uses this knowledge to calculate the coefficients

Briefly, the function you refer to goes something like this:

  • The transfer functions for analog Butterworth filters of any order are already known analytically. It selects the normalized transfer function polynomial for the order you select.

  • The normalized polynomial is scaled appropriately according to the desired cutoff frequency, giving an $s$-domain transfer function that has the response that you want.

  • The $s$-domain prototype filter is then mapped to a discrete-time approximation using the bilinear transform. This is the most commonly-used method of mapping a continuous-time filter to a discrete-time filter.

  • The result of the bilinear transform is a transfer function in the $z$-domain, from which the filter coefficients are easily found by inspection.

$\endgroup$
0
$\begingroup$

For a manually calculation of the coefficients of the IIR-Filter, the Jave speech toolkit has a open source implementation of it.

See Butterworth.java for calculating the coefficients and IIRFilter.java how they are applied on a signal.

Here are additional functions to calculate frequency response, amplitude response and phase response of an arbitary IIR-Filter:

public double amplitudeResponse(final double omega) {
    return ComplexNumber.magnitude(frequencyResponse(omega));
}

public double phaseResponse(final double omega) {
    return ComplexNumber.toEulerPhi(frequencyResponse(omega));
}

private double[] frequencyResponse(final double omega) {
    // z = exp(-j*omega);
    final double[] z = new double[2];
    z[0] = Math.cos(omega);
    z[1] = Math.sin(omega);

    final double[] h = transferFunction(z);
    return h;
}

private @NonNull double[] transferFunction(final double @NonNull [] z) {
    final double[] tmp = new double[2];
    final double[] zp = new double[2];

    final double[] numerator = new double[2];
    numerator[0] = 0;
    numerator[1] = 0;
    for (int k = 0; k < b.length; k++) {
        ComplexNumber.power(z, -k, zp);
        ComplexNumber.multiply(b[k], zp, tmp);
        ComplexNumber.add(numerator, tmp, numerator);
    }

    final double[] denominator = new double[2];
    denominator[0] = 0;
    denominator[1] = 0;
    for (int k = 0; k < a.length; k++) {
        ComplexNumber.power(z, -k, zp);
        ComplexNumber.multiply(a[k], zp, tmp);
        ComplexNumber.add(denominator, tmp, denominator);
    }

    final double[] resesult = new double[2];
    ComplexNumber.divide(numerator, denominator, resesult);
    return resesult;
}

And here are additional unit tests to make it clear how the coefficients $a$ and $b$ are used.

@Test
public void testFilterWithOrderTwoFilter() {
    final float[] x = new float[] { 2f, -3f, 0f, 2f, 0f, 0f, 0f, 0f };
    final double[] b = new double[] { 0.6, 0.4 };
    final double[] a = new double[] { 1.0, 0.5, 0.8 };
    final float[] expectedY = new float[] { 1.2f, -1.6f, -1.36f, 3.16f,
            0.308f, -2.682f, 1.0946f, 1.5983f };

    final float[] actualY = new float[8];
    final IIRFilter filter = new IIRFilter(b, a);
    filter.filter(x, actualY);
    assertArrayEquals(expectedY, actualY, COMPARISON_DELTA);
}

And these unit tests I've used to test against the octave signal package:

@Test
public void testCalcNormalizedFrequency() {
    assertEquals(0.0, calcNormalizedFrequency(0, 44100), COMPARISON_DELTA);
    assertEquals(0.5, calcNormalizedFrequency(44100 / 4, 44100),
            COMPARISON_DELTA);
    // Nyquist
    assertEquals(1.0, calcNormalizedFrequency(44100 / 2, 44100),
            COMPARISON_DELTA);
}

@Test
public void testCreateLowPass() {
    // order 3, cutoff 0.2
    //@formatter:off
    final double[] bExpected = new double[] {
            0.018099, 0.054297, 0.054297, 0.018099 };
    final double[] aExpected = new double[] {
            1.00000, -1.76004, 1.18289, -0.27806 };
    //@formatter:off
    final IIRFilter filter = createLowPass(3, 0.2);
    final double[] bActual = filter.getFeedforwardCoefficients(false);
    final double[] aActual = filter.getFeedbackCoefficients(false);
    assertArrayEquals(bExpected, bActual, COMPARISON_DELTA);
    assertArrayEquals(aExpected, aActual, COMPARISON_DELTA);
}

@Test
public void testCreateHighPass() {
    // order 3, cutoff 0.2
    //@formatter:off
    final double[] bExpected = new double[] {
            0.52762, -1.58287, 1.58287, -0.52762 };
    final double[] aExpected = new double[] {
            1.00000, -1.76004, 1.18289, -0.27806 };
    //@formatter:off
    final IIRFilter filter = createHighPass(3, 0.2);
    final double[] bActual = filter.getFeedforwardCoefficients(false);
    final double[] aActual = filter.getFeedbackCoefficients(false);
    assertArrayEquals(bExpected, bActual, COMPARISON_DELTA);
    assertArrayEquals(aExpected, aActual, COMPARISON_DELTA);
}

@Test
public void testCreateBandPass() {
    // order 3, cutoffLow 0.2, cutoffHigh 0.7
    //@formatter:off
    final double[] bExpected = new double[] {
            0.16667, 0.00000, -0.50000, 0.00000, 0.50000, 0.00000, -0.16667};
    final double[] aExpected = new double[] {
            1.0000e+000, -6.6370e-001, 1.6314e-001, -1.6192e-001,
            3.6596e-001, -7.3744e-002, 1.6374e-017 };
    //@formatter:off
    final IIRFilter filter = createBandPass(3, 0.2, 0.7);
    final double[] bActual = filter.getFeedforwardCoefficients(false);
    final double[] aActual = filter.getFeedbackCoefficients(false);
    assertArrayEquals(bExpected, bActual, COMPARISON_DELTA);
    assertArrayEquals(aExpected, aActual, COMPARISON_DELTA);
}

@Test
public void testCreateBandReject() {
    // order 3, cutoffLow 0.2, cutoffHigh 0.7
    //@formatter:off
    final double[] bExpected = new double[] {
            0.16667, -0.22123, 0.59789, -0.45690, 0.59789, -0.22123, 0.16667 };
    final double[] aExpected = new double[] {
            1.0000e+000, -6.6370e-001, 1.6314e-001, -1.6192e-001,
            3.6596e-001, -7.3744e-002, 1.6374e-017 };
    //@formatter:off
    final IIRFilter filter = createBandReject(3, 0.2, 0.7);
    final double[] bActual = filter.getFeedforwardCoefficients(false);
    final double[] aActual = filter.getFeedbackCoefficients(false);
    assertArrayEquals(bExpected, bActual, COMPARISON_DELTA);
    assertArrayEquals(aExpected, aActual, COMPARISON_DELTA);
}
$\endgroup$
  • $\begingroup$ lotta (potentially very useful) code. less theory. but the code does answer the OP's wish to calculate "manually". $\endgroup$ – robert bristow-johnson Apr 3 '15 at 0:39

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.