Tasks #1253

TTDelay cubic interpolation decision

Added by Nathan Wolek almost 7 years ago. Updated almost 7 years ago.

Status:ClosedStart date:2012-11-27
Priority:NormalDue date:
Assignee:Tim Place% Done:

100%

Category:-Spent time:2.00 hours
Target version:-
Branch:master

Description

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:
http://crca.ucsd.edu/~msp/techniques/latest/book-html/node31.html

I am not sure where the one in TTInterpolate is from. Perhaps Tim could illuminate us?

History

#2 Updated by Trond Lossius almost 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 almost 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.

#4 Updated by Trond Lossius almost 7 years ago

  • Status changed from New to Closed
  • % Done changed from 0 to 100

Applied in changeset commit:b889d1a2484e1048aa24eac7b42d75f7461281d2.

#5 Updated by Trond Lossius almost 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.

Foundation/Library/TTInterpolate:

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.

DSP/Library

This interpolation is based on the following four assumptions:

  1. (f(-1) = x0
  2. f(0) = x1
  3. f(1) = x2
  4. 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).

Conclusion

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 almost 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

Also available in: Atom PDF