[SciPy-user] scipy.signal.chebwin

Ryan May rmay@ou....
Fri Feb 15 15:09:50 CST 2008


Ryan May wrote:
> Kumar Appaiah wrote:
>> On Sat, Feb 09, 2008 at 12:16:04PM +0530, Kumar Appaiah wrote:
>>>> array([-0.16010146, -0.16010146, -0.16010146, -0.16010146, -0.16010146,
>>>>        -0.16010147, -0.16010148, -0.16010149, -0.1601015 , -0.1601015 ,
>>>>        -0.16010145, -0.16010096, -0.16009716, -0.16007336, -0.15994973,
>>>>        -0.15941238, -0.15743963, -0.15127378, -0.13476733, -0.09676449,
>>>>        -0.02138783,  0.10725105,  0.29505955,  0.52638443,  0.7591664 ,
>>>>         0.93452305,  1.        ,  0.93452305,  0.7591664 ,  0.52638443,
>>>>         0.29505955,  0.10725105, -0.02138783, -0.09676449, -0.13476733,
>>>>        -0.15127378, -0.15743963, -0.15941238, -0.15994973, -0.16007336,
>>>>        -0.16009716, -0.16010096, -0.16010145, -0.1601015 , -0.1601015 ,
>>>>        -0.16010149, -0.16010148, -0.16010147, -0.16010146, -0.16010146,
>>>>        -0.16010146, -0.16010146, -0.16010146])
>>>>
>>>> Clearly, all of those negative values are *not* correct.  (And the
>>>> problems are not limited to the numbers above.)  Any ideas?
>>> Let me try to figure it out. Then I'll let you know.
>> I am unable to figure out where the problem could be, though I guess
>> it would have to do with the Chebyshev polynomial evaluation. I could
>> really do with a little help in debugging the chebwin fix. :-)
>>
> Got it!
> 
> It seems scipy.special.chebyt doesn't quite do what we want it to do.
> If I replace the call to chebyt with this function:
> 
> def myT(order, x):
>     retval = N.cosh(order*N.arccosh(N.abs(x)))
>     retval[N.abs(x)<=1] = N.cos(order*N.arccos(x[N.abs(x)<=1]))
>     return retval

This actually needs to be:

def myT(order, x):
    retval = N.zeros_like(x)
    retval[x > 1] = N.cosh(order*N.arccosh(x[x>1]))
    retval[x < -1] = N.cosh(order*N.arccosh(-x[x<-1]))*((-1)*(order%2))
    retval[N.abs(x)<=1] = N.cos(order*N.arccos(x[N.abs(x)<=1]))
    return retval

I missed a problem with odd ordered Tn.  See
http://en.wikipedia.org/wiki/Chebyshev_polynomials.

> I get the correct window values.  I had checked plots of chebyt and the
> analytic definition, and they looked the same, but the actual values
> differed.  For most values they're close, but not identical.  Anybody
> know the difference between the Chebyshev polynomial defined as above
> and the one in scipy.special.chebyt?

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.

Ryan

-- 
Ryan May
Graduate Research Assistant
School of Meteorology
University of Oklahoma


More information about the SciPy-user mailing list