Spline interpolation

[madah]

Hi!

I've been investigating spline interpolation and here is some of my results.
I wanted to find out if there was any way to bridge the gap between the commonly used 4-point cubic spline and 512-point sinc interpolation, both in terms of performance and quality.

All code below is placed in the public domain.


Spline polynomials of third order (aka cubic splines):

• Continuous first derivate everywhere.
• Continuous second derivate on all segments except the last.

Code: Select all

template<class T, class BufT>
T interpolate_spline3_6pt(const T t, const BufT* x)
{
    return (t*(t*(t*(T(17)*(x[2]-x[3]) + T(9)*(x[4]-x[1]) + T(2)*(x[0]-x[5])) +
            (T(20)*(x[1]+x[3]) + T(2)*x[5] - T(4)*x[0] - T(7)*x[4] - T(31)*x[2])) +
            (T(11)*(x[3]-x[1]) + T(2)*(x[0]-x[4])))) / T(14) + x[2];
}

Code: Select all

template<class T, class BufT>
T interpolate_spline3_8pt(const T t, const BufT* x)
{
    return (t*(t*(t*(T(91)*(x[3]-x[4]) + T(45)*(x[5]-x[2]) + T(13)*(x[1]-x[6]) + T(3)*(x[7]-x[0])) +
            (T(106)*(x[2]+x[4]) + T(10)*x[6] + T(6)*x[0] - T(3)*x[7] - T(29)*(x[1]+x[5]) - T(167)*x[3])) +
            (T(61)*(x[4]-x[2]) + T(16)*(x[1]-x[5]) + T(3)*(x[6]-x[0])))) / T(76) + x[3];
}

Spline polynomials of fifth order (aka quintic splines):

• Continuous first, second and third derivate everywhere.
• Continuous fourth derivate on all segments except the last.

Code: Select all

template<class T, class BufT>
T interpolate_spline5_6pt(const T t, const BufT* x)
{
    return (t*(t*(t*(t*(t*(T(841)*(x[3]-x[2]) + T(677)*(x[1]-x[4]) + T(238)*(x[5]-x[0])) +
            (T(2379)*x[2] + T(1559)*x[4] + T(452)*x[0] - T(738)*x[5] - T(1826)*(x[1]+x[3]))) +
            (T(606)*(x[1]-x[3]) + T(500)*(x[5]-x[2]) + T(72)*(x[0]-x[4]))) +
            (T(1820)*(x[1]+x[3]) - T(548)*(x[0]+x[4]) - T(2544)*x[2])) +
            (T(1277)*(x[3]-x[1]) + T(262)*(x[0]-x[4])))) / T(1506) + x[2];
}

Code: Select all

template<class T, class BufT>
T interpolate_spline5_8pt(const T t, const BufT* x)
{
    return (t*(t*(t*(t*(t*(T(6489)*(x[4]-x[3]) + T(4831)*(x[2]-x[5]) + T(3089)*(x[6]-x[1]) + T(1063)*(x[0]-x[7])) +
            (T(17352)*x[3] + T(9062)*(x[1]+x[5]) + T(3204)*x[7] - T(2111)*x[0] - T(6383)*x[6] - T(15093)*(x[2]+x[4]))) +
            (T(6659)*(x[2]-x[4]) + T(5403)*(x[5]-x[1]) + T(2141)*(x[3]-x[7]) + T(45)*(x[6]-x[0]))) +
            (T(17169)*(x[2]+x[4]) + T(2171)*(x[0]+x[6]) - T(5102)*(x[1]+x[5]) - T(28476)*x[3])) +
            (T(13566)*(x[4]-x[2]) + T(4532)*(x[1]-x[5]) + T(1078)*(x[6]-x[0])))) / T(15472) + x[3];
}

Why not have second (or fourth) derivate continuous everywhere?
I had a constraint that there must be unity gain for all values of t, and give straight lines when interpolating stair-steps. This is important since otherwise you'd get "wobbly" lines, which results in a high pitched signal artifact.
This constraint is why the splines can't have their second (or fourth) derivate continuous everywhere, as there is no solution to the equation systems that fulfill both constraints.


Plot of impulse response, together with the 5-point spline from musicdsp:


Close-up of the passband:

[Aleksey Vaneev]

You do not need 512-point sinc interpolation, it's an old concept which is a bit wrong today.

All you really need is oversampling to alleviate the problems of interpolation formulas near the Nyquist frequency. Oversampling shifts these problems beyond area of interest.

I've proved in https://code.google.com/p/r8brain-free-src/ that you just need 38-point sinc interpolator at most given oversampling was done first. The good thing is that oversampling can be performed via fast convolution while there is no way to perform fast fractional interpolation exists.

I would suggest you to start with at least 16-point interpolation, for good results.

[Aleksey Vaneev]

I think if you do 32x oversampling first you can even use plain simple linear interpolation with good results.

[madah]

I think if you do 32x oversampling first you can even use plain simple linear interpolation with good results.

Yes I agree with that. For a softsynth, linear interpolation is often good enough given that the wavetables are large enough (and FFT upsampling of smaller tables is fast and can be done as a preprocessing step during load).

But with 32x oversampling, you either have to store the upsampled parts in an intermediate buffer, or upsample the whole file in one go (requires 32x memory usage). This might be unacceptable in a sampler with a lot of different samples loaded and playing at the same time.

The main benefit with using spline interpolation is that it doesn't require any extra memory storage, or any setup when changing samplerate (e.g. pitch-bend, vibrato, FM, etc). Also you can start playing from anywhere in the sample with 0-sample delay.

But of course, it all depends on your use case. For streaming audio and resampling on-the-fly, using r8brain is probably the best and fastest approach.

[mystran]

But with 32x oversampling, you either have to store the upsampled parts in an intermediate buffer, or upsample the whole file in one go (requires 32x memory usage). This might be unacceptable in a sampler with a lot of different samples loaded and playing at the same time.

Nah, just use polyphase filters and calculate just the samples you need.

I don't even bother with linear interpolation. I just oversample a lot more.

[Aleksey Vaneev]

madah, soft-synths can also use streamed resampling approach, they only need to feed the resampler with enough data (looped if needed). Of course, 32x oversampling is overkill, but 2x or 4x resampling can be used without problems - e.g. by using polyphase filters. Then it's much easier to attain a high quality interpolation with a function of less complexity.

In the future update to r8brain-free-src one can use just 16 taps if needed - this provides 96 dB signal-to-noise ratio (I had to add several fine-tuned windowing functions).

[Aleksey Vaneev]

madah, do you have by any chance formulas for 2nd order interpolation, based on 6 and 8 points?

[madah]

madah, do you have by any chance formulas for 2nd order interpolation, based on 6 and 8 points?

Here are two variants for 6 points, but they aren't very optimal quality-wise:

Code: Select all

template<class T, class BufT>
T spline2_6pt_variant1(const T t, const BufT* x)
{
    return t*(t*(T(2)*(x[1]+x[4]-x[0]-x[5]) - (x[2]+x[3])) +
            (T(2)*(x[0]+x[3]+x[5]-x[1]-x[4]))) + x[2];
}

Code: Select all

template<class T, class BufT>
T spline2_6pt_variant2(const T t, const BufT* x)
{
    return t*(t*((x[1]+x[4]-x[0]-x[2]-x[3]-x[5])) +
            (T(2)*x[3] + (x[0]+x[5]-x[1]-x[4]))) + x[2];
}

As I'm using Maxima for solving the equation systems, I can only get exact solutions.
With only second order, there's too few degrees of freedom to make the curve both smooth and resemble a sinc.

In this case it would probably have been better to use a genetic algorithm. This is something I've been planning to do later.

[Aleksey Vaneev]

Yep, these are not very good.

[FLWrd]

Being interested in this topic as well, I came across Lanczos resampling, which is a variant on sinc. Does anyone have experience with its performance in downsampling?

[madah]

Being interested in this topic as well, I came across Lanczos resampling, which is a variant on sinc. Does anyone have experience with its performance in downsampling?

Lanczos and any sinc-based interpolation will have superior stopband attenuation. The biggest drawback is performance, it's really really slow compared to using a spline. This can be mitigated by using a table-based lookup.

Also, if the sinc is of too small width (too few points), there will be non-unity gain, which might be a problem.

Impulse-response:


Passband:



For already upsampled data, I recommend using the musicdsp "5-point" spline. It has near-zero passband ripple, and the large areas between the "bulbs" are very good at attenuating the aliasing:

[FLWrd]

Thanks.

Wrt Lanczos performance, if you're up or down-sampling by a factor 2 or 4, it should be possible to optimize that by hand by treating more than one sample inside the loop. General resampling seems pretty expensive, indeed.

Thanks for the plots!