[SciPy-Dev] scipy.signal.butter returns complex coefficients

[Charles R Harris]

On Thu, Feb 21, 2013 at 8:48 AM, Pierre Haessig <pierre.haessig at crans.org>wrote:

>  Hi Eric,
>
> Thanks for your feedback !
>
> Le 21/02/2013 00:58, Eric Moore a écrit :
>
> I'm not sure that the correct fix would be, perhaps signal.buttap should
> calculate the poles by calculating half of them, and then conjugating.
> Since I'd bet that the ultimate cause is that we aren't getting exact
> conjugates from the complex exponential.
>
>  That's right. I just made a quick test inspired by signal.buttap code to
> check for this "non conjugate roots issue"
> (
> https://github.com/scipy/scipy/blob/master/scipy/signal/filter_design.py#L1376
> )
>
> N = 4
> n = numpy.arange(1, N + 1)
> p = numpy.exp(1j * (2 * n - 1) / (2.0 * N) * pi) * 1j
> print(p[0])
> print(p[3].conj())
> print(p[0] == p[3].conj())
> print(p[1])
> print(p[2].conj())
> print(p[1] == p[2].conj())
>
> gives :
> (-0.382683432365+0.923879532511j)
> (-0.382683432365+0.923879532511j)
> False
> (-0.923879532511+0.382683432365j)
> (-0.923879532511+0.382683432365j)
> False
>
> Roots are indeed conjugate at about 1e-16 but the equality check in
> np.poly is an exact check.
>
> Now, it may be possible to write signal.buttap in a bit more complex
> fashion to generate true complex conjugate roots but I'm not sure it's a
> good to rely on np.poly to make some "smart conjugate detection". In the
> new polynomial package I'm not sure this is even implemented (
> https://github.com/numpy/numpy/blob/master/numpy/polynomial/polynomial.py#L127).
> From my unexerced eye, the code of
> numpy.polynomial.polynomial.polyfromroots and numpy.poly looks very
> different
>

The multiplications in numpy/polynomial/polynomial are done as a sort of
binary tree, multiplying successively higher degree polynomials together
while trying to keep the roots in each widely separated so as to minimize
numerical error, but it make no attempt to detect conjugate roots. It
probably needs an `asreal` attribute so the user, who knows more about the
problem than the computer, can specify the type of the result.

<snip>

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20130221/a53604a2/attachment.html>