[SciPy-user] scipy.signal.chebwin

[Pauli Virtanen]

Fri, 15 Feb 2008 15:09:50 -0600, Ryan May wrote:
> Ryan May wrote:
[clip]
> In talking with __pv on #scipy, it seems that the problem is likely just
> due to accumulated round off error in using the recurrence relations to
> implement the polynomial.  Since the N-point window function needs an
> (N-1)th degree polynomial, you can see how this begins to really show
> the error.  Is there any reason not to just implement chebyt as above
> within scipy.special?  It would seem to me that a formula would always
> give superior results, but I could (admittedly) be naive about this.

I think the special.chebyt function in scipy does not use recurrence 
relations for evaluating the polynomials at some point (which probably 
would be numerically stable), but instead creates a numpy.poly1d object 
containing the coefficients of the polynomial and uses them to calculate 
results. As the magnitude of the coefficients increases as polynomial 
degree increases, error from cancellations in floating-point arithmetic 
becomes of the order of 1 for n > 60.

This limitation would be useful to mention in the docstrings of functions 
in orthopoly.py

A possible fix could be to implement a way to evaluate polynomials stably 
in numpy.poly1d. (Horner scheme? Using the roots which are anyway known? 
I guess there are several options.)

-- 
Pauli Virtanen