0
$\begingroup$

I am comparing sos (cascaded 2nd order IIR) coefficients from my own implementation of a textbook method versus the coefficients that octave generates.

Here is an snip of the octave code:

fs = 39e3;
fs2 = fs / 2;
f0 = 15000;
N = 8;

[z, p, k] = butter ( N, f0 / fs2, "high" );
octave_sos = zp2sos(z, p, k);

And here is my sos versus octave sos.

>> disp (my_sos);

   0.111341  -0.222681   0.111341   1.000000   1.325538   0.770901
   0.091891  -0.183782   0.091891   1.000000   1.093985   0.461549
   0.081054  -0.162108   0.081054   1.000000   0.964970   0.289186
   0.076191  -0.152383   0.076191   1.000000   0.907076   0.211841
>> disp (octave_sos);

   6.3184e-05  -1.2637e-04   6.3184e-05   1.0000e+00   1.3255e+00   7.7090e-01
   1.0000e+00  -2.0000e+00   1.0000e+00   1.0000e+00   1.0940e+00   4.6155e-01
   1.0000e+00  -2.0000e+00   1.0000e+00   1.0000e+00   9.6497e-01   2.8919e-01
   1.0000e+00  -2.0000e+00   1.0000e+00   1.0000e+00   9.0708e-01   2.1184e-01

I've already verified that the two sets of coefficients have equivalent impulse and frequency response.

My question is what is the simplification that octave is doing on the coefficients? How can some of the coeffecients get transformed into +1, -2? It would be nice to be able to perform the same simplification in order to optimize the filter processing.

$\endgroup$

1 Answer 1

2
$\begingroup$

That's simply a function of gain staging. For simplicity we will ignore real poles and zeros and assume that each section has a complex pole pair ($p_k$) and a complex zero pair ($z_k$). The transfer function can be written as

$$H_k(z) = g_k\frac{(1-z_k\cdot z^{-1})(1-z_k^*\cdot z^{-1})}{(1-p_k\cdot z^{-1})(1-p_k^*\cdot z^{-1})} = \\ g_k\frac{1-2\Re(z_k)\cdot z^{-1} + |z_k|^2 \cdot z^{-2}}{1-2\Re(p_k)\cdot z^{-1} + |p_k|^2 \cdot z^{-2}} \tag{1} = \\ g_k \frac{1+b_{1,k}z^{-1} + b_{2,k}z^{-2}}{1+a_{1,k}z^{-1} + a_{2,k}z^{-2}}$$

where $g_k$ is the gain of the stage. Since the sections are cascaded, the total gain is simply the product of all individual gains, i.e. $g = \prod_k g_k$. So you can factor out the total gain and set the individual section gains to $1$. In practice that means you modify the numerator coefficients as follows

$$g_k = b_{0,k}, b^{'}_{0,k} = 1, b^{'}_{1,k} = b_{1,k}/b_{0,k}, b^{'}_{2,k} = b_{2,k}/b_{0,k} $$

All $b_0$ coefficients become 1. Since a Butterworth filter has all zeros at $z=1$ or $z=-1$, the normalized coefficients become indeed $[1,2,1]$ or $[1,-2,1]$ for all sections.

You can see this by doing

>> [sos1,g]= zp2sos(z, p, k)
sos1 =
    1.0000   -2.0000    1.0000    1.0000    0.9071    0.2118
    1.0000   -2.0000    1.0000    1.0000    0.9650    0.2892
    1.0000   -2.0000    1.0000    1.0000    1.0940    0.4615
    1.0000   -2.0000    1.0000    1.0000    1.3255    0.7709
g =
   6.3184e-05

All section gains are 1 and the overall gain has been factored out as g. The gain can be computed as the product of the first column

>> prod(mysos(:,1))
ans =
   6.3184e-05
>> prod(octave_sos(:,1))
ans =
   6.3184e-05

In order to avoid carrying around an extra gain number, octave uses the simplification of wrapping the entire gain into the first stage, i.e.

>> octave_sos(1,1)
ans =
   6.3184e-05

Note that in most floating point implementations, the individual gain stages are quite flexible, but it is VERY important to manage them properly for fixed point code. This being said, I have certainly come across cases where section gains, section ordering and pole/zero pairing made a big difference in terms of numerical performance.

$\endgroup$
1
  • $\begingroup$ Thanks for the in depth explanation! $\endgroup$
    – mike919
    Commented Apr 28 at 17:33

Your Answer

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

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