Search code examples
pythonscipyinterpolationbounding-boxspline

How do I limit the interpolation region in the InterpolatedUnivariateSpline in Python when given non-uniform samples?


I'm trying to get a nice upsampler using Python when I have non-uniform spaced inputs. Any suggestions would be helpful. I've tried a number of interp functions. Here's an example:

from scipy.interpolate import InterpolatedUnivariateSpline
from numpy import linspace, arange, append
from matplotlib.pyplot import plot
F=[0, 1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,22050]
M=[0.,2.85,2.49,1.65,1.55,1.81,1.35,1.00,1.13,1.58,1.21,0.]
ff=linspace(F[0],F[1],10)
for i in arange(2, len(F)):
    ff=append(ff,linspace(F[i-1],F[i], 10))
aa=InterpolatedUnivariateSpline(x=F,y=M,k=2);
mm=aa(ff)
plot(F,M,'r-o'); plot(ff,mm,'bo'); show()

This is the plot I get:

enter image description here

I need to get interpolated values that don't go below 0. Note that the blue dots go below zero. The red line represents the original F vs. M data. If I use k=1 (piece-wise linear interp) then I get good values as shown here:

aa=InterpolatedUnivariateSpline(x=F,y=M,k=1)
mm=aa(ff); plot(F,M,'r-o');plot(ff,mm,'bo'); show()

enter image description here

The problem is that I need to have a "smooth" interpolation and not the piece-wise value. Does anyone know if the bbox argument in InterpolatedUnivarientSpline helps to fix that? I cant find any documentation on what bbox does. Is there another easier way to accomplish this?

Thanks in advance for any help.


Solution

  • Positivity-preserving interpolation is hard (if it wasn't, there wouldn't be a bunch of papers written about it). The splines of low degree (2, 3) usually do pretty well in this regard, but your data has that large gap in it, and it happens to be at the end of data range, making things worse.

    One solution is to do interpolation in two steps: first upsample the data by piecewise linear interpolation, then interpolate new data with a smooth spline (I'll use cubic spline below, though quadratic also works).

    The gap_size array records how large each gap is, relative to the smallest one. In subsequent loop, uniformly spaced points are replaced in large gaps (those that are at least twice the size of smallest one). The result is F_new, a nearly-uniform better grid that still includes the original points. The corresponding M values for it are generated by a piecewise linear spline. Subsequent cubic interpolation produces a smooth curve that stays positive.

    F = [0, 1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,22050]
    M = [0.,2.85,2.49,1.65,1.55,1.81,1.35,1.00,1.13,1.58,1.21,0.]
    
    gap_size = np.diff(F) // np.diff(F).min()
    F_new = []
    for i in range(len(F)-1):
        F_new.extend(np.linspace(F[i], F[i+1], gap_size[i], endpoint=False))
    F_new.append(F[-1])
    
    pl_spline = InterpolatedUnivariateSpline(F, M, k=1);
    M_new = pl_spline(F_new)
    smooth_spline = InterpolatedUnivariateSpline(F_new, M_new, k=3)
    
    ff = np.linspace(F[0], F[-1], 100)
    plt.plot(F, M, 'ro')
    plt.plot(ff, smooth_spline(ff), 'b')
    plt.show()
    

    spline

    Of course, no tricks can hide the truth that we don't know what happens between 5500 and 22050 (Hz, I presume), the nearly-linear part is just a placeholder.