scipy bspline complex interpolation

bspline_change.py

from numpy import *
from scipy.signal import cubic

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, with N=len(cj)

    Edges are handled using mirror-symmetric boundary conditions.

    """
    newx = (asarray(newx) - x0) / float(dx)
    res = zeros_like(newx, 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
    result = zeros_like(newx, 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

test_cplx_interp.py

#!/usr/bin/env python

from scipy.signal import cspline1d
#from scipy.signal import cspline1d_eval
from bspline_change import cspline1d_eval

from numpy import *

from matplotlib.pylab import *

# create some smoothly varying complex signal to try and interpolate
x=arange(100)
y=zeros(x.shape,dtype=complex64)
T=10.0
f=1.0/T
y=sin(2.0*pi*f*x)+1.0J*cos(2.0*pi*f*x)

# get the cspline transform
cy=cspline1d(y)

# determine some xvalues to interpolate
xnew=arange(1.0,10.0,0.001)
ynew=cspline1d_eval(cy,xnew)
# This warning is generated using current scipy baseline:
#/usr/lib/python2.7/dist-packages/scipy/signal/bsplines.py:348: ComplexWarning:Casting complex values to real discards the imaginary part
#  result += cj[indj] * cubic(newx - thisj)



# throw up some figures comparing these
figure()
plot(xnew,ynew.imag)
figure()
plot(xnew,ynew.real)

figure()
plot(x[1:11],y[1:11].imag)
figure()
plot(x[1:11],y[1:11].real)

show()