[Scipy-tickets] [SciPy] #1260: scipy.signal.cspline1d_eval cannot return complex results

SciPy Trac scipy-tickets@scipy....
Wed Aug 11 09:59:00 CDT 2010

#1260: scipy.signal.cspline1d_eval cannot return complex results
 Reporter:  anielsen      |       Owner:  somebody
     Type:  defect        |      Status:  new     
 Priority:  normal        |   Milestone:  0.9.0   
Component:  scipy.signal  |     Version:  0.7.0   
 Keywords:                |  
 I apologize in advance if I misunderstand what's going on here. I'm new to
 python and scipy.

 I've found that cspline1d_eval cannot return complex results. The forward
 transform (cspline1d) will create complex coefficients when given a
 complex series input. I did some simple tests that confirm these changes
 produce acceptable results for appropriately smooth complex input.

 the code is in scipy/signal/bsplines.py

 I made a simple hack which produces the behavior I'd anticipated by
  res = zeros_like(newx)
  res = zeros(newx.shape,dtype=cj.dtype)
 in two locations, that I marked in the code snippet below.

 I've attached a very simple script which demonstrates the results of using
 these changes to interpolate a simple complex signal.

 def cspline1d_eval(cj, newx, dx=1.0, x0=0):
     """Evaluate a spline at the new set of points.
     dx is the old sample-spacing while x0 was the old origin.

     In other-words the old-sample points (knot-points) for which the cj
     represent spline coefficients were at equally-spaced points of

     oldx = x0 + j*dx  j=0...N-1


     edges are handled using mirror-symmetric boundary conditions.
     newx = (asarray(newx)-x0)/float(dx)
 ### change made here
     #res = zeros_like(newx)
     res = zeros(newx.shape,dtype=cj.dtype)
     if (res.size == 0):
         return res
     N = len(cj)
     cond1 = newx < 0
     cond2 = newx > (N-1)
     cond3 = ~(cond1 | cond2)
     # handle general mirror-symmetry
     res[cond1] = cspline1d_eval(cj, -newx[cond1])
     res[cond2] = cspline1d_eval(cj, 2*(N-1)-newx[cond2])
     newx = newx[cond3]
     if newx.size == 0:
         return res
 ### change made here
     #result = zeros_like(newx)
     result = zeros(newx.shape,dtype=cj.dtype)
     jlower = floor(newx-2).astype(int)+1
     for i in range(4):
         thisj = jlower + i
         indj = thisj.clip(0,N-1) # handle edge cases
         result += cj[indj] * cubic(newx - thisj)
     res[cond3] = result
     return res

Ticket URL: <http://projects.scipy.org/scipy/ticket/1260>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.

More information about the Scipy-tickets mailing list