[SciPy-user] how to make a sinc(x) function - "divide by zero" !

Scott Ransom sransom at nrao.edu
Mon Dec 12 12:04:16 CST 2005


Here is a simple version that I use.  It uses the Taylor 
expansion when the argument is near zero.  Note that this 
is using Numeric and umath right now.  This should be 
trivial to change to new scipy_core (which I am not using yet):

def sinc(xs):
    """
    sinc(xs):
        Return the sinc function [i.e. sin(pi * xs)/(pi * xs)]
            for the values xs.
    """
    pxs = umath.pi*xs
    return Numeric.where(umath.fabs(pxs)<1e-3, 
                         1.0-pxs*pxs/6.0, umath.sin(pxs)/pxs)

Scott


On Monday 12 December 2005 12:39 pm, Arnd Baecker wrote:
> On Mon, 12 Dec 2005, Sebastian Haase wrote:
> > Hi,
> > I was trying to implement a "sinc" [sin(x)/x] in numarray. But the
> > "where"-semantics makes it choke on the x=0: 0/0 case.
> > (This should of course work for x being an array - so "if" is no
> > option) The best I could come up with was:
> >
> > def sinc(r):
> >     na.Error.pushMode(all="ignore")
> >     a = na.where(r, na.divide(na.sin(r),r), 1)
> >     na.Error.popMode()
> >     return a
> >
> > but I still seem to get  a warning...
> >
> > >>> F.sinc(0)
> >
> > Warning: Encountered invalid numeric result(s)  in divide
> > 1.0
> >
> > What is a better way of doing this ?
>
> Yuo could try to use `vectorize` (warning: untested code)
>
> def sinc(x):
>     if x==0.0:                # presumably better to check for small
> x return 0.0            # here ...
>     else:
>        return sin(x)/x
>
> sinc_vectorized=scipy.vectorize(sinc)
>
> HTH,
>
> Arnd
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user

-- 
Scott M. Ransom            Address:  NRAO
Phone:  (434) 296-0320               520 Edgemont Rd.
email:  sransom at nrao.edu             Charlottesville, VA 22903 USA
GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989



More information about the SciPy-user mailing list