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 y_{i}, i = 1...n, where n is the number of samples. Each sample i is associated with an x location (typically, time) given by x_{i} = x_{1} + (i  1) dx, where dx is the sample period, so that the realvalued sample number associated with a given time x is
If the resulting s is an integer number, the y value must be y_{s}. 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:
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:
s_{l} ≡ floor (s) ; s_{r} ≡ s_{l} + 1 
y(s) ≈ y_{l} + (s  s_{l}) · (y_{r}  y_{l}) 
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 y′_{l} and y′_{r} 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} ≈ (y_{r}  y_{l1}) / 2 ; y′_{r} ≈ (y_{r+1}  y_{l}) / 2 
Now that we know y_{l}, y_{r}, y′_{l}, and y′_{r}, we can fit these values with a thirddegree (cubic) polynomial:
As_{l}^{3} + Bs_{l}^{2} + Cs_{l} + D = y_{l} 
As_{r}^{3} + Bs_{r}^{2} + Cs_{r} + D = y_{r} 
3As_{l}^{2} + 2Bs_{l} + C = y′_{l} 
3As_{r}^{2} + 2Bs_{r} + C = y′_{r} 
If we shift the x axis to the left sample point, we can reduce the four equations to
so that the interpolated value y(s) at any point s between s_{l} and s_{r} is estimated as
(y′_{r} + y′_{l}  2y_{r} + 2y_{l}) φ_{l}^{3} + (3y_{r}  3y_{l}  2y′_{l}  y′_{r}) φ_{l}^{2} + y′_{l} φ_{l} + y_{l} 
where φ_{l} ≡ s  s_{l}. Some rearrangement gives
y(s) ≈ y_{l} φ_{r} + y_{r} φ_{l} + 
 φ_{l} φ_{r} [1/2 (y′_{r}  y′_{l}) + (φ_{l}  1/2) (y′_{l} + y′_{r}  2(y_{r}  y_{l}))] 
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 y′_{l} equals y′_{r} and both equal y_{r}  y_{l}, the higherdegree terms are zero.

3. If y′_{l} + y′_{r} equals 2(y_{r}  y_{l}), the thirddegree 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:
s_{l} ≡ floor (s); s_{r} ≡ s_{l} + 1 
φ_{l} ≡ s  s_{l}; φ_{r} ≡ 1  φ_{l} 
y(s) ≈ ∑_{i=1...N} y_{ri} sinc (π(φ_{l}+i1)) (1/2 + 1/2 cos (π(φ_{l}+i1)/(φ_{l}+N))) + 
+ ∑_{i=1...N} y_{l+i} sinc (π(φ_{r}+i1)) (1/2 + 1/2 cos (π(φ_{r}+i1)/(φ_{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 lowpass filtered at the Nyquist frequency. See:

Sound: Get value at time...
Links to this page
© ppgb, September 12, 2020