TTDelay cubic interpolation decision
|Assignee:||Tim Place||% Done:|
|Category:||-||Spent time:||2.00 hours|
I went to replace the cubic interpolation algorithm in TTDelay with the one from TTInterpolate and they seem to be completely different algorithms, but BOTH are simply called cubic. Switching between them shows that they do vary in their results. A web search shows a variety of techniques labelled as cubic. We basically just need to decide: which one do we use?
The one in TTDelay seems to implement this:
I am not sure where the one in TTInterpolate is from. Perhaps Tim could illuminate us?
#1 Updated by Nathan Wolek over 7 years ago
#2 Updated by Trond Lossius over 7 years ago
I'm currently looking into the interpolation algorithms. I don't have any final conclusions yet, but so far the algorithm in DSP seems trustworthy while I'm less convinced about the one in Foundation. But I'm in the process of implementing unit tests for TTInterpolate in Foundation.
#3 Updated by Trond Lossius over 7 years ago
- Branch set to master
I have just added some unit tests for TTInterpolate in Foundation Library (commit f4088a5a5), and they indicate that something is wrong with the cubic interpolator there, at least if it's supposed to work the same way as the interpolator in DSP Library.
#5 Updated by Trond Lossius over 7 years ago
After several days of investigations I have come to the following conclusions:
With the commits done today we now have two differing polynimical interpolation algorithms mplemented. They are now both correctly implemented, but differ in their assumptions.
The assumptions are now thoroughly documented in the header file:
@details This interpolation algorithms calculate the coefficients a, b, c, d of the 3rd order polynomial f(delta) = a*aDelta^3 + b*aDelta^2 + c*aDelta + d = ( (a*aDelta + b )*aDelta + c)*aDelta + d) so that the function fulfill the following four conditions: -# f(0) = x1 @n -# f(1) = x2 @n -# f'(0) = (x2-x0)/2 @n -# f'(1) = (x3-x1)/2 The two last conditions use a symetric estimate of the difference at the end points of the region to interpolate over: 0 ≤ aDelta ≤ 1 These asumptions ensure that the resulting interpolated function, when moving over several subsequent sections, is: -# Continuous (no sudden jump) -# Has a continuous derivative (no break pints with hard edges) However, the 2nd order derivate will generally be discontinuous on the points connecting two sections.
This interpolation is based on the following four assumptions:
- (f(-1) = x0
- f(0) = x1
- f(1) = x2
- f(2) = x3
This ensures that the interpolated function pass through all four points. However, this function has to be expected to have a discontionuous derivative at the connection points (this might appear as hard edges in the resulting interpolated curve).
For audio interpolation purposes the function in Foundation is probably preferable, as it is expected to introduce less artificial spectral noise.
#6 Updated by Trond Lossius over 7 years ago
Here are the calculations involved, documented for posterity:
f(t) = a t^3 + b t^2 + c t + d f'(t) = 3a t^2 + b t + c ---- f(0) = x => d = x (1) --- f'(0) = (y-w)/2 = c => c = -0.5 w + 0.5 y (2) --- f(1) = a + b + c + d = y a + b = y - (-0.5w + 0.5y) - x a + b = y + 0.5w -0.5 y - x 2a + 2b = w - 2x + y (3) --- f'(1) = 3a + 2b + c = (z-x)/2 6a + 4b = z - x - 2c 6a + 4b = z - x - 2 (-0.5 w + 0.5 y) 6a + 4b = z - x + w - y 6a + 4b = w - x - y + z (4) --- Combining (3) and (4) 2a = -w + 3x -3y + z (5) a = -0.5w + 1.5x - 1.5y + 0.5z Combining (3) and (5) 2b = w - 2x + y - ( -w + 3x -3y + z ) 2b = w - 2x + y + w - 3x + 3y - z 2b = 2w - 5x + 4y - z b = w - 2.5x + 2y - 0.5z --- Checking that conditions are satisfied: f(1) = a + b + c + d = ( -0.5w + 1.5x - 1.5y + 0.5z ) + ( w - 2.5x + 2y - 0.5z ) + ( -0.5 w + 0.5 y ) + x = -0.5w + 1.5x - 1.5y + 0.5z + w - 2.5x + 2y - 0.5z - 0.5 w + 0.5 y + x = y f'(1) = 3a + 2b + c = 3( -0.5w + 1.5x - 1.5y + 0.5z ) + 2( w - 2.5x + 2y - 0.5z ) + ( -0.5 w + 0.5 y ) = -1.5w + 4.5x - 4.5y + 1.5z + 2w - 5x + 4y - z - 0.5 w + 0.5 y = 4.5x + 1.5z - 5x - z = -0.5x + 0.5z = (z-x)/2