5
$\begingroup$

The primary source lists Butterworth polynomials in s-domain and provides a link to bilinear transform for digital implementation. But, who needs analog specification in our digital world? Why should everybody get involved into the quite laborious bilinear transform, which, additionally requires frequency pre-warping? Meantime, looking at example bilinear transform in Butterworth design, I see that sampling frequency, T, cancels out, leaving you with a nice polynomial, after the monstrous transform. Why nobody can just conduct those monstrous computations only once and publish the z-polynomials, likewise s-polynomials are provided? Where can I find a table of z-polynomials not to be bogged down with the terrible pre-warping technicalities?

$\endgroup$
3
$\begingroup$

You can't find a table of $z$-domain coefficients that is analogous to the tables of Laplace-domain Butterworth transfer functions because it's impossible. For a second-order Butterworth filter, for example, there is no single corresponding $z$-domain transfer function, in contrast to the continuous-time case.

The reason: discrete-time transfer functions don't cover an infinitely-long frequency axis like a continuous-time system does. When creating a digital approximation to an analog prototype, techniques like the bilinear transform are used to map the frequency response from the $s$-domain transfer function evaluated along the (infinitely long) $j\omega$ axis to the $z$-domain's (finite-length) unit circle. This obviously can't be done exactly, so the various methods for mapping an analog filter to a digital one are approximations.

One factor that plays into the quality of the approximation is the sample rate of the target digital filter. In the paper you linked, you noted that the sampling interval $T$ cancelled out in the calculations of the digital filter coefficients, which is true. However, the sample interval is still required in order to perform the (necessary) frequency prewarping. Recall the prewarping formula:

$$ \omega_p = \frac{2}{T}\tan{\left(\frac{\omega T}{2}\right)} $$

While the rest of the $T$'s may cancel out after performing the substitution $s = \frac{2}{T}\frac{z-1}{z+1}$, the $T$ in the argument to the tangent function cannot go away. Therefore, the resulting coefficients will always be a function of the sample rate $T$, so your hope of creating a simple table of digital Butterworth filter coefficients invariant of sample rate is stymied.

A couple other notes regarding this:

  • The prewarping step is essential to ensure that particular frequencies of interest (i.e. the cutoff frequency for a Butterworth filter) are preserved through the approximation process. Without prewarping, then the stretching of the frequency response that occurs when mapping from an infinitely-long frequency axis to the unit circle would cause the cutoff to appear at a different frequency than desired in the digital filter.

  • The cookbook linked in niaren's answer also demonstrates the dependence upon the sample rate: you have to trudge through the equations a bit, but you'll see factors in the listed coefficients like $\cos(\omega_0)$ and so on (this is really just a way of slightly hiding the prewarping process). The cookbook defines $\omega_0$ as: $$ \omega_0 = \frac{2\pi f_0}{F_s} $$ so when it all comes down to it, you'll have to define your sample rate and prewarp any frequencies of interest in order to design a Butterworth filter approximation. It's not particularly difficult, as the cookbook shows; I'm not sure what the motivation for avoiding those steps is.

$\endgroup$
  • $\begingroup$ 1. Don't you consider the coefficients for e.g. the lowpass filter in the cookbook as coefficients for a Butterworth filter? 2. Do you seriously mean that you don't see the motivation behind making the cookbook? $\endgroup$ – niaren Jun 1 '13 at 15:28
  • 1
    $\begingroup$ @niaren: They are coefficients for a Butterworth filter, but they are dependent upon $\omega_0$, which is dependent upon the sample rate $F_s$. I didn't say that there's no reason for the cookbook; I was merely responding to the OP's question on why there isn't a single table of Butterworth $z$-domain transfer functions like there is for the $s$ domain. In a discrete-time filter, the coefficients will always be dependent upon the filter's sample rate, so there isn't a canonical set of numerical coefficients for digital Butterworth filters like there is for the continuous-time case. $\endgroup$ – Jason R Jun 1 '13 at 16:31
  • $\begingroup$ I was going to write something but found out I had to think some more and accidentally I pressed enter. Found out I couldn't revert. Anyway, nice post to think about. $\endgroup$ – niaren Jun 1 '13 at 17:40
  • $\begingroup$ Given a sample, then there is no table for Butterworth z-domain transfer functions because that table will not be very useful. The reason being that the mapping from coefficients in the normalized z-transfer function to the desired z-transfer function is not trivial. If it had been just a linear mapping of the frequency variable I think tables would be around even though they would depend on sample rate. So maybe it is not so much the actual dependence upon Fs as it is the convoluted mapping from coeffient in normalized to coefficient in denormalized... $\endgroup$ – niaren Jun 2 '13 at 19:41
2
$\begingroup$

This has been done for biquad designs: cookbook. Is this what you're looking for?

Edit: Let's try and make a z-transfer function table for Butterworth. We have to choose some normalized cutoff frequency $\omega_n $(in radians pr second). Let us choose e.g. $0.1\pi$.

N=2: $\frac{0.02008 + 0.04017z^{-1} + 0.02008z^{-2}}{1 - 1.561z^{-1} + 0.6414z^{-2}}$

N=3: $\frac{0.002898 + 0.008695z^{-1} + 0.008695z^{-2}+0.002898z^{-3}}{1 - 2.374z^{-1} + 1.929z^{-2} -0.5321z^{-3}}$

N=4: $\frac{0.0004166 + 0.001666 ^{-1} + 0.002500 ^{-2}+0.001666 ^{-3}+0.0004166^{-4}}{1 -3.181z^{-1} + 3.861z^{-2} -2.112z^{-3} + 0.4383z^{-4}}$

In order for this table to be useful we should be able to use it to design filters for any desired cutoff frequency. If we could linearly map the frequency variable e.g. $p = cz$, as in the analog case the table would be just as useful as the analog table. However, the mapping to be useful must be one that maps the unit circle in the complex $z$ plane onto itself. We choose the simplest mapping

$z^{-1} \rightarrow \frac{-\alpha + z^{-1}}{1 - \alpha z^{-1}} $ (1)

To obtain the coefficients of the denormalized z-transfer function we have to compute $\alpha$ in the following way

$\alpha = \frac{\sin\left(\frac{\omega_n- \omega_c}{2}\right)}{\sin\left(\frac{\omega_n+ \omega_c}{2}\right)}$

where $\omega_c$ is the new cutoff frequency of the Butterworth LPF. In order to compute the coefficients in the denormalized z-transfer function (1) must be used in the normalized z-transfer function. If we do this for the case N=2, we obtain

$\bar{b}_0 = b_0-\alpha b_1+b_2\alpha^2 \\ \bar{b}_1 = -2\alpha b_0+b_1+\alpha^2b_1-2\alpha b_2 \\ \bar{b}_2 = b_0\alpha^2-\alpha b_1+b_2 \\ \bar{a}_0 = 1-\alpha a_1+a_2 \alpha^2 \\ \bar{a}_1 = a_1+a_1\alpha^2 -2\alpha a_2-2\alpha \\ \bar{a}_2 = -\alpha a_1+a_2+\alpha^2$

Not only is it painful to derive these expressions but the amount of computations that must be done in a design implementation seem to be not that much simpler than starting out from the analog normalized case.

Edit2: Here is a 'proof' showing that the z-transfer function table approach gives the same coefficients as the bilinear transformation from an analog prototype.

Bilinear transformation to any normalized angular frequency, $\omega_c$,

$s \rightarrow \frac{1}{\tan(\omega_c/2)}\frac{1-z^{-1}}{1+z^{-1}} \quad (1)$

In the z-transfer function table approach the bilinar transformation is done for one particular frequency, $\omega_n$

$s \rightarrow \frac{1}{\tan(\omega_n/2)}\frac{1-z^{-1}}{1+z^{-1}} \quad (2) $

Then frequency warping is performed to the desired cutoff frequency by frequency warping

$z^{-1} \rightarrow \frac{-\alpha+z^{-1}}{1-\alpha z^{-1}} \quad (3)$

Inserting (3) into (2) gives

$s \rightarrow \frac{1}{\tan(\omega_n/2)}\frac{1+\alpha}{1-\alpha}\frac{1-z^{-1}}{1+z^{-1}}$

Substitution of $\alpha = \frac{\sin(\omega_c - \omega_n)}{\sin(\omega_c + \omega_n)}$ into the above equation gives (1).

I will have to look at the complexity of the two methods later.

$\endgroup$
  • 1
    $\begingroup$ I have computed the the presampling koefficient, $tg=\tan(\pi f_{cutoff}/f_{sample})$, and succeeded turning $H(z) = {tg^2 \over ({1-z^{-1} \over 1+z^{-1} })^2 +\sqrt 2 tg ({1-z^{-1} \over 1+z^{-1} }) + tg^2}$ into $H(z) = {tg^2(1+2z^{-1}+z^{-2}) \over 1+\sqrt 2 tg + tg^2 + 2(tg^2-1)z^{-1} + (1-\sqrt 2 tg+ tg^2)z^{-2}}$, as www.robots.ox.ac.uk/~sjrob/Teaching/SP/l6.pdf teaches us. It seems that I need to compute 3 sums to compute every coefficient. Your result, though different, is identical. Yet, it seems much simpler than making the laborious conversion that I undergone, every time. $\endgroup$ – Val Jun 5 '13 at 19:21
  • $\begingroup$ I have updated the post. Maybe that helps a little in the understanding. Regarding the complexity, this is definitely something I'm going to take a look at as soon as possible. $\endgroup$ – niaren Jun 7 '13 at 13:14
0
$\begingroup$

Just a function + factors computed by me for 4th order low-pass butterworth filter. C++. Immensely big thank you to @Val above who showed how to compute 2nd order filter.

yb and xb in the example are just storage in order it would be possible to make processing by separate buffers. freq_hz is a cut-off frequency in hertz.

    static void butterworth4(float * const __restrict signal, int const signal_len, float (& __restrict yb)[4], float (& __restrict xb)[4],
                             float const freq_hz, float const sample_rate)
    {
            float const f = tan(M_PI * freq_hz / sample_rate);      //prewarped frequency
            float const f2 = f * f;
            float const f3 = f2 * f;
            float const f4 = f2 * f2;
            float const m = -2.f * cosf(5.f * M_PI / 8.f);
            float const n = -2.f * cosf(7.f * M_PI / 8.f);

            float const a0 = 1.f + (m + n) * f + (2.f + n * m) * f2 + (m + n) * f3 + f4;
            float const a1 = (-4.f - 2.f * (n + m) * f + 2.f * (m + n) * f3 + 4.f * f4) / a0;
            float const a2 = (6.f - 2.f * (2 + m * n) * f2 + 6.f * f4) / a0;
            float const a3 = (-4.f + 2.f * (m + n) * f - 2.f * (m + n) * f3 + 4.f * f4) / a0;
            float const a4 = (1.f - (n + m) * f + (2.f + m * n) * f2 - (m + n) * f3 + f4) / a0;

            for (int i = 0; i < signal_len; ++i) {
                    float const num = f4 / a0 * (signal[i] + 4.f*xb[0] + 6.f*xb[1] + 4.f*xb[2] + xb[3]);
                    float const den = -a1 * yb[0] - a2 * yb[1] - a3 * yb[2] - a4 * yb[3];
                    float const y = num + den;
                    yb[3] = yb[2]; yb[2] = yb[1]; yb[1] = yb[0]; yb[0] = y;
                    xb[3] = xb[2]; xb[2] = xb[1]; xb[1] = xb[0]; xb[0] = signal[i];
                    signal[i] = y;
            }
    }

    static void butterworth2(float * const __restrict signal, int const signal_len, float (& __restrict yb)[2], float (& __restrict xb)[2],
                             float const freq_hz, float const sample_rate)
    {
            float const f = tan(M_PI * freq_hz/sample_rate);
            float const f2 = f * f;

            float const sq2 = sqrtf(2.f);
            float const a0 = 1 + sq2 * f + f2;
            float const a1 = 2.f * (f2 - 1.f) / a0;
            float const a2 = (1 - sq2 * f + f2) / a0;

            for (int i = 0; i < signal_len; ++i) {
                    float const num = f2 / a0 * (signal[i] + 2.f*xb[0] + xb[1]);
                    float const den = -a1 * yb[0] - a2 * yb[1];
                    float const y = num + den;
                    yb[1] = yb[0]; yb[0] = y;
                    xb[1] = xb[0]; xb[0] = signal[i];
                    signal[i] = y;
            }
    }
$\endgroup$

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.