CubicSplineFast.java (7249B)
1 /* 2 * _______ _____ _____ _____ 3 * |__ __| | __ \ / ____| __ \ 4 * | | __ _ _ __ ___ ___ ___| | | | (___ | |__) | 5 * | |/ _` | '__/ __|/ _ \/ __| | | |\___ \| ___/ 6 * | | (_| | | \__ \ (_) \__ \ |__| |____) | | 7 * |_|\__,_|_| |___/\___/|___/_____/|_____/|_| 8 * 9 * ------------------------------------------------------------- 10 * 11 * TarsosDSP is developed by Joren Six at IPEM, University Ghent 12 * 13 * ------------------------------------------------------------- 14 * 15 * Info: http://0110.be/tag/TarsosDSP 16 * Github: https://github.com/JorenSix/TarsosDSP 17 * Releases: http://0110.be/releases/TarsosDSP/ 18 * 19 * TarsosDSP includes modified source code by various authors, 20 * for credits and info, see README. 21 * 22 */ 23 24 25 /********************************************************** 26 * 27 * Class CubicSplineFast 28 * 29 * Class for performing an interpolation using a cubic spline 30 * setTabulatedArrays and interpolate adapted, with modification to 31 * an object-oriented approach, from Numerical Recipes in C (http://www.nr.com/) 32 * Stripped down version of CubicSpline - all data checks have been removed for faster running 33 * 34 * 35 * WRITTEN BY: Dr Michael Thomas Flanagan 36 * 37 * DATE: 26 December 2009 (Stripped down version of CubicSpline: May 2002 - 31 October 2009) 38 * UPDATE: 14 January 2010 39 * 40 * DOCUMENTATION: 41 * See Michael Thomas Flanagan's Java library on-LineWavelet web page: 42 * http://www.ee.ucl.ac.uk/~mflanaga/java/CubicSplineFast.html 43 * http://www.ee.ucl.ac.uk/~mflanaga/java/ 44 * 45 * Copyright (c) 2002 - 2010 Michael Thomas Flanagan 46 * 47 * PERMISSION TO COPY: 48 * 49 * Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee, 50 * provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies 51 * and associated documentation or publications. 52 * 53 * Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice, 54 * this list of conditions and the following disclaimer and requires written permission from the Michael Thomas Flanagan: 55 * 56 * Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and 57 * the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission 58 * from the Michael Thomas Flanagan: 59 * 60 * Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose. 61 * Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software 62 * or its derivatives. 63 * 64 ***************************************************************************************/ 65 66 67 package be.tarsos.dsp.util; 68 69 public class CubicSplineFast{ 70 71 private int nPoints = 0; // no. of tabulated points 72 private double[] y = null; // y=f(x) tabulated function 73 private double[] x = null; // x in tabulated function f(x) 74 private double[] d2ydx2 = null; // second derivatives of y 75 76 // Constructors 77 // Constructor with data arrays initialised to arrays x and y 78 public CubicSplineFast(double[] x, double[] y){ 79 this.nPoints=x.length; 80 this.x = new double[nPoints]; 81 this.y = new double[nPoints]; 82 this.d2ydx2 = new double[nPoints]; 83 for(int i=0; i<this.nPoints; i++){ 84 this.x[i]=x[i]; 85 this.y[i]=y[i]; 86 } 87 this.calcDeriv(); 88 } 89 90 // Constructor with data arrays initialised to zero 91 // Primarily for use by BiCubicSplineFast 92 public CubicSplineFast(int nPoints){ 93 this.nPoints=nPoints; 94 this.x = new double[nPoints]; 95 this.y = new double[nPoints]; 96 this.d2ydx2 = new double[nPoints]; 97 } 98 99 // METHODS 100 101 // Resets the x y data arrays - primarily for use in BiCubicSplineFast 102 public void resetData(double[] x, double[] y){ 103 for(int i=0; i<this.nPoints; i++){ 104 this.x[i]=x[i]; 105 this.y[i]=y[i]; 106 } 107 } 108 109 // Returns a new CubicSplineFast setting array lengths to n and all array values to zero with natural spline default 110 // Primarily for use in BiCubicSplineFast 111 public static CubicSplineFast zero(int n){ 112 if(n<3)throw new IllegalArgumentException("A minimum of three data points is needed"); 113 CubicSplineFast aa = new CubicSplineFast(n); 114 return aa; 115 } 116 117 // Create a one dimensional array of cubic spline objects of length n each of array length m 118 // Primarily for use in BiCubicSplineFast 119 public static CubicSplineFast[] oneDarray(int n, int m){ 120 CubicSplineFast[] a =new CubicSplineFast[n]; 121 for(int i=0; i<n; i++){ 122 a[i]=CubicSplineFast.zero(m); 123 } 124 return a; 125 } 126 127 128 // Calculates the second derivatives of the tabulated function 129 // for use by the cubic spline interpolation method (.interpolate) 130 // This method follows the procedure in Numerical Methods C language procedure for calculating second derivatives 131 public void calcDeriv(){ 132 double p=0.0D,qn=0.0D,sig=0.0D,un=0.0D; 133 double[] u = new double[nPoints]; 134 135 d2ydx2[0]=u[0]=0.0; 136 for(int i=1;i<=this.nPoints-2;i++){ 137 sig=(this.x[i]-this.x[i-1])/(this.x[i+1]-this.x[i-1]); 138 p=sig*this.d2ydx2[i-1]+2.0; 139 this.d2ydx2[i]=(sig-1.0)/p; 140 u[i]=(this.y[i+1]-this.y[i])/(this.x[i+1]-this.x[i]) - (this.y[i]-this.y[i-1])/(this.x[i]-this.x[i-1]); 141 u[i]=(6.0*u[i]/(this.x[i+1]-this.x[i-1])-sig*u[i-1])/p; 142 } 143 144 qn=un=0.0; 145 this.d2ydx2[this.nPoints-1]=(un-qn*u[this.nPoints-2])/(qn*this.d2ydx2[this.nPoints-2]+1.0); 146 for(int k=this.nPoints-2;k>=0;k--){ 147 this.d2ydx2[k]=this.d2ydx2[k]*this.d2ydx2[k+1]+u[k]; 148 } 149 } 150 151 // INTERPOLATE 152 // Returns an interpolated value of y for a value of x from a tabulated function y=f(x) 153 // after the data has been entered via a constructor. 154 // The derivatives are calculated, bt calcDeriv(), on the first call to this method ands are 155 // then stored for use on all subsequent calls 156 public double interpolate(double xx){ 157 158 double h=0.0D,b=0.0D,a=0.0D, yy=0.0D; 159 int k=0; 160 int klo=0; 161 int khi=this.nPoints-1; 162 while (khi-klo > 1){ 163 k=(khi+klo) >> 1; 164 if(this.x[k] > xx){ 165 khi=k; 166 } 167 else{ 168 klo=k; 169 } 170 } 171 h=this.x[khi]-this.x[klo]; 172 173 if (h == 0.0){ 174 throw new IllegalArgumentException("Two values of x are identical: point "+klo+ " ("+this.x[klo]+") and point "+khi+ " ("+this.x[khi]+")" ); 175 } 176 else{ 177 a=(this.x[khi]-xx)/h; 178 b=(xx-this.x[klo])/h; 179 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; 180 } 181 return yy; 182 } 183 }