Cherry Audio Forums

Code: Select all

import java.util.Arrays;

// cubic spline interpolation
 
public class CubicSpline
{
   // constructor for immutable instance
   // precomputes tangents for subsequent calls to interpolate()
   // yInputs is an array of at least two Y values that represent the control points (aka nodes)
   // the x values are implicit integers equidistant from 0 to the length of yInputs - 1
   // if monotone is true the tangents are adjusted to prevent non-monotonic artifacts
   // i.e. overshoot that doesn't seem appropriate given the input data
   // but in practice it can be safely set to false to speed things up
   // note the yInput data does not need to be monotonic
   
	public CubicSpline( double[] yInputs, boolean monotone )
	{
	   // n is the number of control points
		int n = yInputs.length;
		
		// make our own copy to decouple from the outside world...
		y = Arrays.copyOf( yInputs, n );

		double[] d = new double[ n - 1 ];      // temp array of slopes
		m = new double[ n ];    // array of tangents

		// calc slopes...
		for( int i = 0; i < n - 1; i++ )
			d[ i ] = y[ i + 1 ] - y[ i ];    // the change in y

		// calc tangets...
		m[ 0 ] = d[ 0 ];
		for( int i = 1; i < n - 1; i++ )
			m[ i ] = ( d[ i - 1 ] + d[ i ] ) * 0.5;   // the average of the slopes
		m[ n - 1 ] = d[ n - 2 ];

      // optionally modify tangents to preserve monotonicity...
      if( monotone )
      {
         for( int i = 0; i < n - 1; i++ )
         {
            if( d[ i ] == 0 )
            {
               // zero slope
               m[ i ] = 0;
               m[ i + 1 ] = 0;
            }
            else
            {
               // non-zero slope
               double a = m[ i ] / d[ i ];
               double b = m[ i + 1 ] / d[ i ];
               double h = Math.hypot( a, b );
               if( h > 9 )
               {
                  // adjust the tangents...
                  double t = 3 / h;
                  m[ i ] = t * a * d[ i ];
                  m[ i + 1 ] = t * b * d[ i ];
               }
            }
		   }
      }
   }


   // interpolate f( x ) for precomputed values using cubic Hermite spline interpolation 
   // 0 <= x <= number of y values - 1
   // out of bounds x handled safely by clamping
	// results pass exactly through every control point
	
	public double interpolate( double x )
	{
		// handle limits...
		if( x <= 0 )
			return y[ 0 ];
			
		int maxIndex = y.length - 1;			
		if( x >= maxIndex )
			return y[ maxIndex ];

      // array index is integer part of x...
      int i = (int) Math.floor( x );

      // difference is fractional part of x...
		double t = x - i;
		
		
		// compute the cubic Hermite spline interpolation...
		
		// h00( t ) * y0 +
		// h10( t ) * m0 +
		// h01( t ) * y1 +
		// h11( t ) * m1
		
		// where the basis functions are the following polynomials...
		// h00( t ) == 2t^3 - 3t^2 + 1
		// h10( t ) == t^3 - 2t^2 + t
		// h01( t ) == -2t^3 + 3t^2
		// h11( t ) == t^3 - t^2
		
		// these can be rearranged as follows...
		// h00( t ) == ( 1 + 2 * t ) * ( 1 - t ) ^ 2
		// h10( t ) == t * ( 1 - t) ^ 2
		// h01( t ) == t ^ 2 * ( 3 - 2 * t )
		// h11( t ) == t ^ 2 * ( t - 1 )
		
		// and coded efficiently like so...
	
		return ( y[ i ] * ( 1 + 2 * t ) + m[ i ] * t ) * ( 1 - t ) * ( 1 - t )
				+ ( y[ i + 1 ] * ( 3 - 2 * t ) + m[ i + 1 ] * ( t - 1 ) ) * t * t;
	}
	
   private double[] m;     // array of tangents
   private double[] y;     // array of y values

}