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)
Improve this question
$\endgroup$
10
  • 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, 2012 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, 2012 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, 2012 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, 2012 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, 2012 at 17:38

2 Answers 2

Reset to default
7
$\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
Improve this answer
$\endgroup$
3
  • $\begingroup$ All right, I will follow your last advice and use numerical stable FIR filters. :) Good and detailed answer though. $\endgroup$
    – alfa
    Apr 8, 2012 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, 2012 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, 2012 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.

Improve this answer
$\endgroup$
3
  • 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, 2012 at 21:54
  • $\begingroup$ Nice. I learnt something new :) $\endgroup$
    – learnvst
    Apr 11, 2012 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, 2014 at 19:57

Your Answer

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

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