5
$\begingroup$

This is a follow-up question of How do I apply a Chebishev filter? I designed a filter in Scipy:

import scipy.signal as signal
fs = 240
signal.cheby1(8, 0.05, [0.1/(fs/2), 10.0/(fs/2)],btype='band', analog=0, output='ba')
/usr/lib/python2.7/dist-packages/scipy/signal/filter_design.py:268: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  "results may be meaningless", BadCoefficients)
(array([ -3.64482771e-08,  -3.98755170e-23,   1.27568970e-07,
     3.98755170e-23,  -2.55137939e-07,  -2.09346464e-22,
     3.18922424e-07,   1.48536301e-21,  -2.55137939e-07,
    -2.79128619e-22,   1.27568970e-07,  -3.98755170e-23,
    -3.64482771e-08,   0.00000000e+00,   4.55603463e-09]),
array([  1.00000000e+00,  -1.53699059e+01,   1.10865293e+02,
    -4.98175837e+02,   1.56087517e+03,  -3.61582779e+03,
     6.40630344e+03,  -8.85525782e+03,   9.65121227e+03,
    -8.32140081e+03,   5.65721533e+03,  -3.00062986e+03,
     1.21729735e+03,  -3.65135628e+02,   7.63720484e+01,
    -9.95191544e+00,   6.08648759e-01]))

The first array should contain the b coefficients, the second array should contain the b coefficients of the filter.

The result of this filter is oscillating and seems to be wrong. Hilmar wrote that this filter "is a numerically challenging filter, since you have poles very very close to the unit circle" and I must "break the filter down into second order sections and apply those sequentially".

Since I have no background in signal processing and found nothing about that anywhere I wanted to ask here how I can do that.

edit: The zero-pole-gain ('zpk') representation is

signal.cheby1(8, 0.05, [0.1/(fs/2), 10.0/(fs/2)],btype='band', analog=0, output='zpk')
/usr/lib/python2.7/dist-packages/scipy/signal/filter_design.py:268: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
"results may be meaningless", BadCoefficients)
(array([ 0.95041122+0.63504447j,  0.95041122-0.63504447j,
    0.77688699+0.32179713j,  0.77688699-0.32179713j,
    0.72152481+0.14352021j,  0.72152481-0.14352021j,
    0.70710678+0.j        , -0.95041122+0.63504447j,
   -0.95041122-0.63504447j, -0.77688699+0.32179713j,
   -0.77688699-0.32179713j, -0.72152481+0.14352021j,
   -0.72152481-0.14352021j, -0.70710678+0.j        ]), array([ 1.19049737+0.j        ,  1.16303520+0.1084667j ,
    1.16303520-0.1084667j ,  1.08874047+0.19530989j,
    1.08874047-0.19530989j,  0.99220878+0.25677213j,
    0.99220878-0.25677213j,  0.93226615+0.2765259j ,
    0.93226615-0.2765259j ,  0.87852354+0.22994371j,
    0.87852354-0.22994371j,  0.83333317+0.15993377j,
    0.83333317-0.15993377j,  0.80434290+0.08130729j,
    0.80434290-0.08130729j,  0.79450811+0.j        ]), -3.6448277062436296e-08)
$\endgroup$
  • 3
    $\begingroup$ Find the roots of the transfer function polynomial (or poles of the rational function) which will be in conjugate pairs $\alpha$ and $\alpha^*$, and then form $(z^{-1}-\alpha)(z^{-1}-\alpha^*)$ to get second-order sections. $\endgroup$ – Dilip Sarwate Apr 8 '12 at 0:05
  • 1
    $\begingroup$ @Dilip: I respectfully disagree here. If you are in polynomial form, it's already too late. In this specific case you cannot recover the poles/zeros, at least not with double precision math). You need to start with poles & zeros $\endgroup$ – Hilmar Apr 8 '12 at 12:55
  • 2
    $\begingroup$ @Hilmar Numerical issues (which are very severe here) aside, are you claiming that it is mathematically impossible to create cascaded second-order sections from a filter specification as a polynomial or a rational function? Or is it that it is far, far better to have designed the filter in terms of poles and zeroes to begin with and avoid the hassle of root finding? If the latter, I am in full agreement. My comment was a partial answer to the OP's question which asked how to find second-order sections given his present state. I agree that it would be better to have started from elsewhere. $\endgroup$ – Dilip Sarwate Apr 8 '12 at 13:15
  • 2
    $\begingroup$ @Dilip: my concern is purely numerical, but very real. The noise gain from a polynomial coefficient to a root location can be arbitrarily high. Even if you have arbitrary precision math available, you may not be able to tell what precision is actually needed. $\endgroup$ – Hilmar Apr 8 '12 at 13:49
  • 2
    $\begingroup$ @alfa: the zpk results are bad too. There poles in there with a magnitude larger than 1, so that's garbage. It looks like you can't use the cheby1() function for that filter. $\endgroup$ – Hilmar Apr 8 '12 at 17:38
6
$\begingroup$

Below the results as sos matrix and gain. Every row in the matrix is a single second order section in the coefficient order of b0, b1, b2, a0, a1, a2. Obviously we need a0=1, and the sections are normalized so that b0=1 as well. The overall gain accumulated is in the single number "gain". The fact that this gain is actually in the order of 10^-9 shows that this is a numerically tricky beast.

Here is how this works.

  1. Design the filter in poles and zeros (zp or zpk), NOT in polynomial form (ba). I imagine that the cheby1() function has an output option for that. If not, you can't use it.
  2. The poles come in conjugate complex pairs, the zeros in this case are real, but you can think of them as a conjugate complex pair with a phase of 0.
  3. Group one pole pair each with one zero pair. Start with the pole pair that has the highest magnitude and find that zero pair that's closest to it. This step is really important!!
  4. Turn the poles and zeros into polynomial in 1/z of the shape (1-p/z)*(1-p'/z). Let's say the complex zero/pole is x+j*y. Then b0/a0 coefficients are always 1, the b1/a1 coefficients are -2*x and the b2/a2 coefficients are x^2+y^2.
  5. Repeat over all poles. Each conjugate pole pair results in one second order section.

A word of caution, if I may: If none of things that I have listed above mean anything to you, than you probably should not be doing this. This is not trivial and many things can go wrong in non-obvious ways. Depending on what you are planning to do with the results, you really want to make sure that your code is fully vetted and properly tested. Unless you have a good idea on how to test this, I would not base anything mission critical on the results.

    sos =
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.824520683146258   0.835012167790563
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.825351261746257   0.854128332434601
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.845749373866514   0.899792317177862
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.890243419701735   0.963689345198141
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.989176997083584   0.989218209927750
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.996496707908275   0.996511870590091
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.998714416931744   0.998722644911976
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.999653030269727   0.999659262704001
gain =
    4.556034632804534e-009
$\endgroup$
  • $\begingroup$ All right, I will follow your last advice and use numerical stable FIR filters. :) Good and detailed answer though. $\endgroup$ – alfa Apr 8 '12 at 17:04
  • $\begingroup$ FIR may not work either. You probably need around 10000 tabs to get this transfer function. This is a very aggressive filter spec. How big is your data vector? $\endgroup$ – Hilmar Apr 8 '12 at 17:35
  • $\begingroup$ I use this low pass filter at the moment: signal.firwin(30+1, 0.1/fs, window='hamming'), no high pass filter. My data vector has actually 7794 (240 Hz) samples (* 64 channels) but I use only 160 (0.667 seconds) samples. At least the signal is smooth now. $\endgroup$ – alfa Apr 8 '12 at 17:41
2
$\begingroup$

If you do not care about how to arrive at the answer, do not need to update the filter coefficients continuously, and have access to Matlab, then you can use the tf2sos command.

[sos,gain] = tf2sos(b,a)

...will give you the output in Hilmar's answer.

$\endgroup$
  • 2
    $\begingroup$ Not in this case. The double precision math that MATLAB is using is not good enough to find the roots properly. [sos, gain] = zp2sos(z,p,k), will do the trick. $\endgroup$ – Hilmar Apr 11 '12 at 21:54
  • $\begingroup$ Nice. I learnt something new :) $\endgroup$ – learnvst Apr 11 '12 at 22:13
  • $\begingroup$ Note that this has been implemented in scipy, too, and will appear in future versions: github.com/scipy/scipy/pull/3717 $\endgroup$ – endolith Jul 7 '14 at 19:57

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.