core.be.tarsos.dsp.util.CubicSplineFast Maven / Gradle / Ivy
The newest version!
/*
* _______ _____ _____ _____
* |__ __| | __ \ / ____| __ \
* | | __ _ _ __ ___ ___ ___| | | | (___ | |__) |
* | |/ _` | '__/ __|/ _ \/ __| | | |\___ \| ___/
* | | (_| | | \__ \ (_) \__ \ |__| |____) | |
* |_|\__,_|_| |___/\___/|___/_____/|_____/|_|
*
* -------------------------------------------------------------
*
* TarsosDSP is developed by Joren Six at IPEM, University Ghent
*
* -------------------------------------------------------------
*
* Info: http://0110.be/tag/TarsosDSP
* Github: https://github.com/JorenSix/TarsosDSP
* Releases: http://0110.be/releases/TarsosDSP/
*
* TarsosDSP includes modified source code by various authors,
* for credits and info, see README.
*
*/
/**********************************************************
*
* Class CubicSplineFast
*
* Class for performing an interpolation using a cubic spline
* setTabulatedArrays and interpolate adapted, with modification to
* an object-oriented approach, from Numerical Recipes in C (http://www.nr.com/)
* Stripped down version of CubicSpline - all data checks have been removed for faster running
*
*
* WRITTEN BY: Dr Michael Thomas Flanagan
*
* DATE: 26 December 2009 (Stripped down version of CubicSpline: May 2002 - 31 October 2009)
* UPDATE: 14 January 2010
*
* DOCUMENTATION:
* See Michael Thomas Flanagan's Java library on-LineWavelet web page:
* http://www.ee.ucl.ac.uk/~mflanaga/java/CubicSplineFast.html
* http://www.ee.ucl.ac.uk/~mflanaga/java/
*
* Copyright (c) 2002 - 2010 Michael Thomas Flanagan
*
* PERMISSION TO COPY:
*
* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
* and associated documentation or publications.
*
* Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice,
* this list of conditions and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
*
* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission
* from the Michael Thomas Flanagan:
*
* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
* or its derivatives.
*
***************************************************************************************/
package be.tarsos.dsp.util;
public class CubicSplineFast{
private int nPoints = 0; // no. of tabulated points
private double[] y = null; // y=f(x) tabulated function
private double[] x = null; // x in tabulated function f(x)
private double[] d2ydx2 = null; // second derivatives of y
// Constructors
// Constructor with data arrays initialised to arrays x and y
public CubicSplineFast(double[] x, double[] y){
this.nPoints=x.length;
this.x = new double[nPoints];
this.y = new double[nPoints];
this.d2ydx2 = new double[nPoints];
for(int i=0; i=0;k--){
this.d2ydx2[k]=this.d2ydx2[k]*this.d2ydx2[k+1]+u[k];
}
}
// INTERPOLATE
// Returns an interpolated value of y for a value of x from a tabulated function y=f(x)
// after the data has been entered via a constructor.
// The derivatives are calculated, bt calcDeriv(), on the first call to this method ands are
// then stored for use on all subsequent calls
public double interpolate(double xx){
double h=0.0D,b=0.0D,a=0.0D, yy=0.0D;
int k=0;
int klo=0;
int khi=this.nPoints-1;
while (khi-klo > 1){
k=(khi+klo) >> 1;
if(this.x[k] > xx){
khi=k;
}
else{
klo=k;
}
}
h=this.x[khi]-this.x[klo];
if (h == 0.0){
throw new IllegalArgumentException("Two values of x are identical: point "+klo+ " ("+this.x[klo]+") and point "+khi+ " ("+this.x[khi]+")" );
}
else{
a=(this.x[khi]-xx)/h;
b=(xx-this.x[klo])/h;
yy=a*this.y[klo]+b*this.y[khi]+((a*a*a-a)*this.d2ydx2[klo]+(b*b*b-b)*this.d2ydx2[khi])*(h*h)/6.0;
}
return yy;
}
}