vector value interpolation

An algorithm for estimating the value of a sampled signal at a specified location.

### Overview

The signal is described with the sequence yi, i = 1...n, where n is the number of samples. Each sample i is associated with an x location (typically, time) given by xi = x1 + (i - 1) dx, where dx is the sample period, so that the real-valued sample number associated with a given time x is

 s = (x - x1) / dx + 1

If the resulting s is an integer number, the y value must be ys. Otherwise, the estimated y value y(s) must be interpolated from nearby values of y. The precision of the result depends on the interpolation method of this algorithm.

### 1. Lowest precision: round to sample

If the interpolation method is “nearest”, we take the value of the nearest point in either direction:

 near ≡ round (s)
 y(s) ≈ ynear

### 2. Middle precision: linear interpolation

If you know or assume that the function that underlies your points is continuous, the "rounding" interpolation would be poor, because the rounded value would abruptly change at the centres between the sample points.

For a linear interpolation, therefore, we use the attested values on both sides (left and right) of s:

 sl ≡ floor (s) ; sr ≡ sl + 1
 y(s) ≈ yl + (s - sl) · (yr - yl)

where floor (x) computes the greatest integer not greater than x. This interpolation is continuous.

### 3. Higher precision: cubic interpolation

If you know or assume that the function that underlies your points is smooth, i.e. its derivative is defined for every x, linear interpolation would probably be poor, because the derivative of the interpolated function would abruptly change at every sample point.

The next higher interpolation (“cubic”), therefore, is differentiable at sample points. To enforce this, we define the derivatives yl and yr at the left and right sample points on the basis of their immediate neighbours (i.e., the algorithm needs four sample points), perhaps by a parabolic interpolation through these three points. A parabolic interpolation has the advantage that the extrema will be computed correctly if the underlying function can be approximated by a parabola near its extremes (see vector peak interpolation).

Because the derivative of a parabolic function is a linear function of x, the derivatives at the left and right sample points are simply estimated as

 y′l ≈ (yr - yl-1) / 2 ; y′r ≈ (yr+1 - yl) / 2

Now that we know yl, yr, yl, and yr, we can fit these values with a third-degree (cubic) polynomial:

 Asl3 + Bsl2 + Csl + D = yl
 Asr3 + Bsr2 + Csr + D = yr
 3Asl2 + 2Bsl + C = y′l
 3Asr2 + 2Bsr + C = y′r

If we shift the x axis to the left sample point, we can reduce the four equations to

 D = yl
 A + B + C + D = yr
 C = y′l
 3A + 2B + C = y′r

so that the interpolated value y(s) at any point s between sl and sr is estimated as

 (y′r + y′l - 2yr + 2yl) φl3 + (3yr - 3yl - 2y′l - y′r) φl2 + y′l φl + yl

where φls - sl. Some rearrangement gives

 y(s) ≈ yl φr + yr φl +
 - φl φr [1/2 (y′r - y′l) + (φl - 1/2) (y′l + y′r - 2(yr - yl))]

where φr ≡ 1 - φl. From this formula we see:

1. The first two terms define the linear interpolation.
2. If the underlying function is linear, so that yl equals yr and both equal yr - yl, the higher-degree terms are zero.
3. If yl + yr equals 2(yr - yl), the third-degree term is zero, so that the interpolated function is parabolic.
4. At the left and right points, one of the φ is 0 and the other is 1, so that at these boundary points, y is computed with exact precision.

### 4. Highest precision: sinc interpolation

If the interpolation method is “sinc70” or “sinc700”, the algorithm assumes that the signal is a sum of sinc functions, so that a number of points (the interpolation depth: 70 or 700) on each side of s must be taken into account.

Because the interpolation depth must be finite, the sum of sinc functions is multiplied by a Hanning window:

 sl ≡ floor (s); sr ≡ sl + 1
 φl ≡ s - sl; φr ≡ 1 - φl
 y(s) ≈ ∑i=1...N yr-i sinc (π(φl+i-1)) (1/2 + 1/2 cos (π(φl+i-1)/(φl+N))) +
 + ∑i=1...N yl+i sinc (π(φr+i-1)) (1/2 + 1/2 cos (π(φr+i-1)/(φr+N)))

where the sinc function is defined as

 sinc (0) ≡ 1; sinc (x) ≡ sin x / x for x ≠ 0

If s is less than the interpolation depth or greater than n + 1 minus the interpolation depth, the depth is reduced accordingly.

This method is appropriate for signals that result from sampling a continuous signal after it has been low-pass filtered at the Nyquist frequency. See:

Sound: Get value at time...