Mike is a senior software engineer at Computerized Medical Systems in St. Louis, Missouri. He can be reached at [email protected].
When collecting or computing discrete data (particularly time-varying data), the resulting resolution often is limited by influences such as computational burden, A/D sampling rates, storage limitations, and transmit rates. This can be a problem when you are interested in the data that describes the minimum and/or maximum values within the continuous function. Examples of this include the determination of the true peak of a computed-signal cross-correlation profile and the reconstruction of the precise peaks and valleys of an environmental measurement.
Determining these minima and maxima, collectively referred to as "extrema," requires interpolation between the known discrete data points. However, if the interpolation resolution is too course, the true extrema might be overlooked; if the resolution is too fine, the computation takes entirely too long. William Press et al. point out a number of algorithms for more efficiently finding the extrema of a data set in Numerical Recipes in C, Second Edition (Cambridge University Press, 1988). These algorithms typically consist of bracketed searches and are often iterative.
In this article, I'll present an alternative approach derived from the well-established equation for cubic-spline interpolation. (For background information on splines, see the accompanying text box entitled "Splinal Tap, or What's My Spline?") This approach performs a direct, noniterative determination of the extrema of a discrete data set.
Algorithm Overview
The cubic-spline extrema algorithm computes the relative extrema of the continuous function that describes the discrete data set. A relative or local extremum is the highest or lowest value within a finite portion of the input data set. A global or absolute extremum is the highest or lowest value within the entire data set. Therefore, given all the relative extrema, the absolute extrema can be determined.
The method uses all the given data to compute second derivatives at each point (also called "knots"). The space between each knot (an "interval") is analyzed in the cubic-spline sense to determine if and where extrema exist. This process directly yields the x values, with no iterative searching. The x values are then used to compute their corresponding y values. This entire process is fast and accurate, and its implementation is quite robust.
Algorithm Derivation
The goal of cubic-spline interpolation is to derive an interpolation formula in which the first and second derivatives of the spline polynomials are equal at the knots. This results in a formula with interval splines that touch at the knots and exhibit a smooth transition from interval to adjacent interval.
Given a data set described by the general function yj = y(xj), the linear interpolation in the interval between xi and xi+1 can be expressed as Figure 1 (a), where " denotes the second derivative and Figure 1 (b) is true. The first derivative of Figure 1 (a) is denoted as dy/dx and is solved as Figure 1 (c).
The derivation of these equations is available in a number of books, including Press et al. However, if you take this further and set Figure 1 (c) to zero, then a new equation can be derived that represents the maximization and minimization of Figure 1 (a). This new equation will allow identification of the points at which y remains constant with respect to finite changes in x.
Finally, Figure 1 (c) can be expressed in the quadratic form ax2 + bx + c = 0, such that x can be solved to yield the cubic-spline extrema equation of Figure 1 (d).
Using these cubic-extrema quadratic coefficients and a quadratic-root solver yields the candidate extrema x1 and x2. If one or both of these abscissa lie within the current interval of examination xi to xi+1, then the candidate is a valid abscissa value at which an ordinate extremum exists. If neither lies within the interval, then no extremum lies within the interval.
Implementation
Listing One is cubic_io.c, a data input and output program that reads and parses an input ASCII data file, calls FindCubicExtrema(), and prints the returned relative extrema. The format of the ASCII data is one x, y point per line, where the x and y are either real or integer and are separated by a comma. Each x value should be greater than the previous one.
In main(), the data file is opened and the number of points is determined so that the input data arrays x_in[] and y_in[] can be allocated to a length of num_pnts. Then the file pointer is reset to the beginning of the file and each line is read and parsed into the x_in and y_in arrays. Also allocated is the structure which will contain the first of a singly linked list of computed extrema. This structure is defined in cubic.h (Listing Two). Each allocation should be checked for success before continuing. However, these checks and their associated error handling (for example, message posting, memory freeing returning an error code) are not shown in the listings.
Because I am not using a doubly linked list that could provide pointers to the previous data, I return to the first structure in the list by saving its address in the pointer variable first_extr prior to calling FindCubicExtrema().
FindCubicExtrema() is then called and passed the number of input data points, the input x and y arrays, and the first output extr structure. If no extremum is found, the function returns a Failure and the user is notified. If the function returns Success, the linked-list pointer is reset to the first one and each computed extremum is listed as the linked list is stepped through until a Null pointer is reached.
All relative extrema are returned--not just the absolute extrema--because no assumptions are made about the user's application. For example, computing the average period of the input data would require all extrema, and examining the rate of decay of an underdamped oscillatory function would require all maxima. If only an absolute extremum were desired (for a cross-correlation peak for example), then it would be very simple to loop through all the relative extrema to find the absolute maximum. In fact, the loop that lists the relative extrema could be easily modified to perform this task. For these reasons, the implementation of FindCubicExtrema() is generalized such that all relative extrema are computed and returned to the calling program.
Cubic-Extrema Computation
The primary routine that computes the data extrema is FindCubicExtrema() (see Listing Three), which performs the following:
- 1. Calls ComputeSecDerivs() to compute the second derivatives of each interval.
- 2. Calculates the cubic extrema quadratic coefficients and calls FindQuadRoots() to compute
the roots of the quadratic equation.
- 3. Determines if the roots, which are candidate extremum abscissa, lie within the current
interval.
- 4. Calls ComputeY() to determine the corresponding ordinate value for each valid extremum
abscissa.
Computing the Second Derivatives
The derivative form of the cubic-spline equation in Figure 1 (c) has everything required to perform the calculation of the quadratic coefficients except the second derivatives of the input data ordinates, denoted as y". Therefore, these must first be computed using the routine ComputeSecDerivs().
This computation is performed as the solution of a system of N spline equations in N unknowns of the general form a11x1 + a12x2 + ... + a1NxN = b1 through aN1x1 + aN2x2 + ... + aNNxN = bN. These equations may be represented in matrix form as Ax = b, where A represents the matrix of as. This system of equations is a symmetric tridiagonal form that is easily solved with a degenerate version of Gaussian elimination. In the software implementation, however, different variables are used to avoid confusion with the quadratic coefficients of FindCubicExtrema(). The matrix variables denoted as a are instead represented in the code as main_diag[] and diag[], which correspond respectively to the main diagonal and off-diagonal elements. The bs are represented as right[] to indicate the right side of the equations. You are solving for the second derivative variables x, which are denoted as sec_deriv[].
The computation of the second derivatives in ComputeSecDerivs() involves division by the diagonals, which are computed from the abscissa data x[]. The requirement that the abscissa values be increasing implies that the data must be unique, which is the case for the problem I have posed. If the abscissa did not change but the ordinate did, then the data would represent an impulse function and the corresponding derivative would be infinity.
If these abscissa data were nonunique, a divide-by-zero error would occur. Therefore, to practice defensive programming, the values of diag[] and main_diag are asserted to ensure that they are not zero. A debug-environment assertion approach is used instead of an "if...then return Failure" approach because this is a nonrecoverable error and, according to the rules I have established, this condition should theoretically never occur. If it does, something has gone awry in allowing this data to get through earlier defenses. If a "ship" version of this routine were not guaranteed unique data, additional error handling would be required beyond the assert().
Computing the Quadratic Roots
Because the second derivatives at the knots are now available, the cubic-extrema quadratic coefficients are easily solved using the equations in Figure 1 (d). These coefficients are then passed to the routine FindQuadRoots() to determine candidate abscissa extrema. Because one or more of the coefficients may be very small, the results are computed using a quadratic solver that avoids associated accuracy errors.
The software implementation of FindQuadRoots() includes checks for conditions that could result in disastrous operations--taking the square root of a negative number or dividing by zero. Either one of these indicates early-on that the interval under scrutiny does not have extrema and thus a Failure is returned. Furthermore, if only one of the divisors is zero, then its respective status (Failure1 or Failure2) will be returned to FindCubicExtrema() so that it knows which of the two roots is a valid candidate. If all these tests pass, then a general status is returned to indicate that both are candidate extrema in the associated interval. Assertions are not used because these conditions can occur in a final product and are not the result of improper programming.
Performing Bounds Checking
If FindQuadRoots() returns a status other than Failure, then at least one of the roots is a valid abscissa candidate. The roots' return status is tested, and if they are accepted for further testing, bounds checking is performed to determine if the root lies within the current interval. If it does, then it is a valid root and thus an abscissa location at which an extremum lies.
Returning the specific Failure1 or Failure2 status from FindQuadRoots() avoids finding an invalid root in the current interval. For instance, if the variable a were computed to be zero and thus x1 were not computed, the previously computed value for x1 would still persist in FindCubicExtrema(). It would then be quite possible for this value to fall outside of the previous interval bounds but now be within the current interval bounds and thus be incorrectly selected as a valid root. Additionally, the variables x1 and x2 could not simply be set to a known value, such as zero, to indicate failure because that may also be a valid abscissa within the next interval.
Extremum-Ordinate Computation
If a valid abscissa has been found, it is then utilized to compute the corresponding ordinate value with ComputeY(). The second derivatives have already been computed so this computation is easy using the standard equation for a cubic spline.
If one or more roots were found in FindCubicExtrema(), Success is returned to main(). Thus, no interrogation of the list is required unless valid data exists. Furthermore, since the first structure was allocated in main(), this is the only way to determine if the structure contained valid data, outside of passing back the number of extrema. If no extrema were found, then the next pointer would never have been set to Null. Therefore, prior to returning Failure in FindCubicExtrema(), the next pointer is Nulled so that the loop in main() responsible for freeing the linked list can recognize the end of the list.
Test Cases and Results
Ensuring that software is accurate and consistent requires that it be validated using input data that yield calculable and thus comparable outputs. These different data sets represent test cases that provide performance verification and "stress" testing. I developed several test cases for testing, debugging, and validating the cubic-spline extrema algorithm and its software implementation.
Four-point Data. The first data set is simple, consisting of the four data points plotted in Figure 2, with a polynomial interpolating curve overlaid for visual aid. Symmetry implies that a maximum is expected at an x value of 0.0 and a y value of approximately 1.1. This data set permits validation that the program can handle negative input data and solve for an extremum at an abscissa of zero. I ran the program using these data and the extremum computed were as expected:
Expected: x = 0.0, y~1.1 Computed: x = 0.00, y = 1.15
Three-point Data. The last three of the four data points from the previous data were then used to illustrate two other important attributes of the algorithm and software implementation. It is intuitively apparent that two data points cannot yield an extremum--only a straight line. At least three are required, therefore representing a minimal data set. This "stress" test also shows that nothing in the algorithm or software precludes an extremum from being located in the first or last interval. This data set was input and yielded the results of 0.0774 and 1.096 for the x and y extremum, respectively. These represent errors of less than 10 percent, which is quite satisfactory given the small amount of data provided to the algorithm.
Trajectory data. The next data set is an extension of the first in that it yields a single maximum. However, this seven-point data set is derived from the equation expressing the trajectory of a projectile. Therefore, the expected abscissa and ordinate values are calculable. The data plotted in Figure 3 was hand-calculated based upon a projectile shot at 1500 feet/second at an initial angle of +45 degrees and a position sampling at 10-second intervals until ground impact.
The standard trajectory equations for the maximum height and corresponding distance were then used to hand calculate the expected extremum for comparison to the algorithm results.
Expected: x = 34,937.89, y = 17,468.94 Computed: x = 34,896.04, y = 17,469.07
The x error was only about 0.10 percent and the y error was negligible.
Cubic splines operate best when provided exact data. But what if the data were approximate due to some real-world measurement errors. To evaluate the impact of this on the algorithm, the trajectory data was rounded to the nearest 20 feet. This provided a good simulation of the errors that might be expected if the global positioning system (GPS) were used to track the projectile with an accuracy of approximately +/-18 feet. This rounded data set was input to the program and yielded the results of x = 34,960.77 and y = 17,485.27, which each translate to errors of less than 0.10 percent.
System response data. The next data set is from control-system theory. Given the step response
of the underdamped second-order system y(t) = 1 - 1.414 e-t cos(4t
- 45
) and evaluated from t = 0 to 3.5 seconds in 0.25 steps, the 15 points are yielded as in Figure 4.
As illustrated, this data set should yield both maxima and minima. An interpolating curve is not shown in Figure 4 because even a single fifth order polynomial is not of high-enough order to characterize the function represented by the data. However, by computing the first derivative of the given function y(t), its expected theoretical extrema can be hand calculated. Five expected extrema were found using this function-dependent calculation. Supplying the data set to the cubic extrema algorithm yielded the following results:
Expected: t = 0.135, y(t) = -1.986 t = 0.920, y(t) = 1.547 t = 1.706, y(t) = 0.751 t = 2.491, y(t) = 1.114 t = 3.277, y(t) = 0.949 Computed: t = none, y(t) = none t = 0.915, y(t) = 1.551 t = 1.707, y(t) = 0.751 t = 2.492, y(t) = 1.114 t = 3.279, y(t) = 0.949
The first extremum was not found, although the cause is obvious from observation of Figure 4. The data supplied to the algorithm indicated no trend toward the first minimum. Thus, the algorithm could not have been expected to locate this extremum without more finely sampled data. However, of the extrema detected, the largest extremum error noted in the time axis was only 0.5 percent and the maximum y(t) error was about 0.25 percent.
Two-Root Data. The last data I'll present provide the validation of several algorithm/software characteristics which represent a stress-test situation. These are the ability to perform with irregularly sampled and/or missing data, to find two extrema within a single interval, and to find very subtle extrema. This combination is a rare occurrence that represents a worst-case scenario. The eight-point data set depicted in Figure 5 comprises samples of the function f(x) = x3 - 2x2 + x + 1.
Note that the 4th interval is quite large compared to the others. Manually taking the derivative of this function shows both a maximum and a minimum in this interval. Intuitively, it makes sense that there can be no more than two extrema in an interval. The combination of irregular sampling and two subtle extrema make this situation challenging and uncommon. The function's expected and resulting data are:
Expected: x = 0.333, f(x) = 1.148 x = 1.000, f(x) = 1.000 Computed: x = 0.318, f(x) = 1.147 x = 1.030, f(x) = 1.007
These data indicate abscissa errors of less than 5 percent and ordinate errors of less than 1 percent. This is very successful given the circumstances. However, it has been shown empirically that if the first data point were not provided to the program, then the first extremum would not be detected. This is because although a root is calculated in the third interval, its value falls within the second interval due to extremum subtlety, and the root is therefore invalidated by the bounds checking. In such situations, the more data, the better.
Miscellaneous
I conducted several additional tests to validate the performance of the algorithm and code. These included a degenerate case of both a positively and a negatively sloped straight ramp. For these cases, no extrema exist and none were detected. I also tested a data set that has decreasing abscissa values, which is a violation of our increasing-data condition. Using this data showed that while no extrema were calculable with this data, the software did not fault, either. Additionally, I used a data set which itself contains an extremum. This was important to ensure that the root computations and bounds-checking logic would not preclude an extremum that coincides with an interval knot.
Conclusions
The cubic-spline extrema algorithm effectively determines the relative extrema of a given data set. It can accommodate negative input data and solve for an extremum at an abscissa of zero. An extremum can be found with just three data points, and nothing in the algorithm or software precludes an extremum from being located in the first or last interval. Two extremum can be detected per interval, subtle extremum are calculable, and accurate results are achievable with and without exact input values. Finally, the higher the data resolution and accuracy, the more accurate the computed extrema will be.
Splinal Tap, or What's My Spline ?
Curve or data fitting finds a function that matches a set of observed values. The primary methods are interpolation and least-squares data fitting. Interpolation assumes the given data to be exact and attempts to find a smooth function that passes through all the data points. Least-squares fitting assumes the given data to be approximate and determines a function that passes through or near all the points as well as possible.
A disadvantage of some interpolators, such as polynomial interpolators, is that a single polynomial function may not accurately satisfy all of the data points, unless of higher order, and thus error may be introduced. Spline interpolation, on the other hand, computes a different polynomial in each interval. These polynomials are computed such that every point, or knot, is exactly intersected and a smooth transition exists from interval to interval.
The term "spline" originates from drafting, where flexible pieces of wood (splines) were used to draw smooth curves by bending them between knots. The shape assumed by the spline between the knots is a third-degree (or cubic) polynomial.
--M.J.C.
Figure 1: Equations for deriving the algorithm.
Figure 2: Simple maximum data.
Figure 3: Trajectory data.
Figure 4: Underdamped second order system.
Figure 5: Two-root function.
Listing One
/* File: cubic_io.c -- Input/output program for the FindCubicExtrema() routine.
Reads an ascii file of x,y data, calls the routine, and lists the
returned extrema. -- by Mike J. Courtney
*/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "cubic.h"
main()
{
/* the following are required for main() */
unsigned int i; /* array indices */
unsigned int num_extr; /* number of extrema found */
BYTE len; /* input data parsing variable */
char data_file[80]; /* ascii data file name */
char str[40], line[80]; /* input data parsing variables */
FILE *fp; /* ascii data file pointer */
struct point *first_extr; /* pointer to the first extreme structure */
/* the following are required for FindCubicExtrema() */
unsigned int num_pnts; /* number of x,y pairs */
float *x_in, *y_in; /* arrays of x and y input data values */
struct point *extr; /* pointer to current extreme structure */
/* open input data file and determine number of input points */
printf("\n Enter name of input file: ");
gets(data_file);
if((fp = fopen(data_file, "r")) != NULL)
{
num_pnts = 0;
while(!feof(fp))
{
fgets(line, 80, fp);
num_pnts++;
}
num_pnts -= 1;
printf("\n The number of input points was %d.",num_pnts);
/* allocate the input data arrays */
x_in = malloc(num_pnts * sizeof(float));
y_in = malloc(num_pnts * sizeof(float));
/* read in the each data line and parse into x and y arrays */
rewind(fp);
for (i = 0; i < num_pnts; i++)
{
/* get a line */
fgets(line, 80, fp);
len = strcspn(line,",");
/* get the x value */
strncpy(str, line, len);
/* NULL terminate the x value */
str[len] = '\0';
x_in[i] = atof(str);
/* get the y value */
strcpy(str, line+len+1);
y_in[i] = atof(str);
}
fclose(fp);
/* allocate first structure of linked list of output extrema */
extr = (struct point *)malloc(sizeof(struct point));
/* save the address of the first structure in the list */
first_extr = extr;
/* call the routine that computes the extrema */
if (FindCubicExtrema(num_pnts, x_in, y_in, extr) == FAILURE)
printf("\n\7 No extrema found !\n");
else
{
/* print the linked list of extrema */
printf("\n\n Relative extrema computed:");
extr = first_extr;
num_extr = 0;
while (extr)
{
printf("\n %d. y = %f at x = %f", num_extr+1, extr->y, extr->x);
extr = extr->next;
num_extr++;
}
}
free(x_in);
free(y_in);
/* free the linked list of extreme structures */
do
{
/* point to first structure */
extr = first_extr;
/* save address of next extreme structure so that when we free the
current one, we still have a pointer to the next one */
first_extr = extr->next;
free(extr);
}
while (first_extr != NULL);
}
else
printf("\n Couldn't open %s !", data_file);
}
Listing Two
/* File: cubic.h. -- #defines, typedefs, and prototypes for cubic_io.c &
cubic_extrema.c. -- by Mike J. Courtney
*/
/* status defines */
#define SUCCESS 1
#define FAILURE 0
#define FAILURE1 -1
#define FAILURE2 -2
/* flag defines */
#define TRUE 1
#define FALSE 0
/* typedefs */
typedef int BOOL;
typedef unsigned char BYTE;
/* structures */
struct point {
struct point *next;
float x;
float y;
};
/* prototypes */
BOOL ComputeSecDerivs (unsigned int, float *, float *, float *);
BOOL FindCubicExtrema (unsigned int, float *,float *, struct point *);
BOOL FindQuadRoots (float, float, float, float *, float *);
void ComputeY (unsigned int, float *, float, float *, float *, float *);
Listing Three
/* File: cubic_extrema.c -- Contains FindCubicExtrema() and supporting routines
that implement the cubic spline extrema algorithm. Given a set of x,y data
points, determine in the cubic spline sense the relative extrema of the
function describing the data. -- by Mike J. Courtney
*/
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include "cubic.h"
/* Primary routine that implements the cubic spline extrema algorithm. Calls
ComputeSecDerivs() to compute the second derivatives, computes quadratic
coefficients, calls FindQuadRoots() to solve quadratic roots, determines if
roots are valid abscissa, and calls ComputeY() to compute ordinates.
*/
BOOL FindCubicExtrema (
unsigned int num_pnts, /* input - address of number of x & y points */
float *x, /* input - array of x values */
float *y, /* input - array of y values */
struct point *extr) /* output - singly linked list of extrema */
{
float a, b, c; /* coefficients of quadratic equation */
float x1, x2; /* roots of quadratic equation to be computed */
unsigned int i; /* array index */
float *sec_deriv; /* array of second derivatives of each data interval */
BOOL root_stat; /* computation status returned by FindRoots() */
BOOL valid_flag; /* TRUE if at least one valid root found */
BOOL first_flag; /* TRUE if current root is the first one */
/* allocate array for second derivatives */
sec_deriv = malloc(num_pnts * sizeof(float));
/* compute the second derivatives */
ComputeSecDerivs(num_pnts, x, y, sec_deriv);
/* initialize extrema flags */
valid_flag = FALSE;
first_flag = TRUE;
/* loop through all the input points and find the extrema */
for (i = 0; i < num_pnts - 1; i++)
{
/* compute the quadratic coefficients */
a = 3.0 * (sec_deriv[i+1] - sec_deriv[i]);
b = 6.0 * (x[i+1] * sec_deriv[i] - x[i] * sec_deriv[i+1]);
c = 6.0 * (y[i+1] - y[i]) - (2.0 * x[i+1] * x[i+1] - x[i] * x[i] + 2 *
x[i] * x[i+1]) * sec_deriv[i];
c -= (x[i+1] * x[i+1] - 2.0 * x[i] * x[i+1] - 2.0 * x[i] * x[i]) *
sec_deriv[i+1];
/* determine the roots of the cubic extrema quadratic equation */
root_stat = FindQuadRoots(a, b, c, &x1, &x2);
if (root_stat != FAILURE)
{
/* if root x1 was calculated successfully */
if (root_stat != FAILURE1)
{
/* Determine if root is within the interval */
if ((x1 > x[i]) && (x1 < x[i+1]))
{
/* first root (extremum) */
if (first_flag == TRUE)
first_flag = FALSE;
/* beyond first valid root so allocate next extremum structure */
else
{
extr->next = (struct point *)malloc(sizeof(struct point));
extr = extr->next;
}
extr->next = NULL;
extr->x = x1;
/* compute the corresponding value of y at the extreme x value */
ComputeY(i, x, extr->x, y, sec_deriv, &extr->y);
valid_flag = TRUE;
}
}
/* if root x2 was calculated successfully */
if (root_stat != FAILURE2)
{
/* Determine if root is within the current interval */
if ((x2 > x[i]) && (x2 < x[i+1]))
{
/* first root (extremum) */
if (first_flag == TRUE)
first_flag = FALSE;
/* beyond first valid root so allocate next extremum structure */
else
{
extr->next = (struct point *)malloc(sizeof(struct point));
extr = extr->next;
}
extr->next = NULL;
extr->x = x2;
/* compute the corresponding value of y at the extreme x value */
ComputeY(i, x, extr->x, y, sec_deriv, &extr->y);
valid_flag = TRUE;
}
}
} /* end of if(root_stat ! = FAILURE) */
} /* end of for(i) */
free(sec_deriv);
if (valid_flag == TRUE)
return SUCCESS;
else
{
/* Set next to NULL just in case it was not set in the loop - this is
so that free loop will operate properly upon return */
extr->next = NULL;
return FAILURE;
}
}
/* Use input x,y data to form tridiagonal matrix and compute second
derivatives of function in the cubic spline sense. */
BOOL ComputeSecDerivs (
unsigned int num_pnts, /* input - number of x & y points */
float *x, /* input - array of x values */
float *y, /* input - array of y values */
float *sec_deriv) /* output - array of 2nd derivatives of intervals */
{
unsigned int i; /* index */
float ftemp; /* temporary float */
float *main_diag; /* ptr to matrix main diagonal array */
float *diag; /* ptr to matrix diagonal array */
float *right; /* ptr to array of right sides of matrix equations */
main_diag = malloc((num_pnts - 2) * sizeof(float));
diag = malloc((num_pnts - 2) * sizeof(float));
right = malloc((num_pnts - 2) * sizeof(float));
/* compute the matrix main and off-diagonal values */
/* even though the calling program is suppose to have guaranteed that the
x values are increasing, assert that neither of the diagonal
differences are zero to avoid a divide by zero condition */
for (i = 1; i < num_pnts - 1; i++)
{
main_diag[i-1] = 2.0 * (x[i+1] - x[i-1]);
assert(main_diag[i-1] > 0);
}
for (i = 0; i < num_pnts - 1; i++)
{
diag[i-1] = x[i+1] - x[i];
assert(diag[i-1] > 0);
}
/* compute right hand side of equation */
for (i = 1; i < num_pnts - 1; i++)
right[i-1] = 6.0 * ((y[i+1]-y[i])/diag[i-1]-(y[i]-y[i-1])/diag[i-2]);
/* forward eliminate tridiagonal */
sec_deriv[0] = 0.0;
sec_deriv[num_pnts - 1] = 0.0;
for (i = 1; i < num_pnts - 2; i++)
{
ftemp = diag[i] / main_diag[i];
right[i] -= (right[i-1] * ftemp);
main_diag[i] -= (diag[i-1] * ftemp);
}
/* backward substitution to solve for second derivative at each knot */
for (i = num_pnts - 2; i > 0; i--)
sec_deriv[i] = (right[i-1] - diag[i-1] * sec_deriv[i+1]) / main_diag[i-1];
free(main_diag);
free(diag);
free(right);
return SUCCESS;
}
/* Solve for roots x1 and x2 of a quadratic equation of the form
a * (x * x) + b * x + c = 0 using the following formula x1 = d / a and
x2 = c / d, where d = -0.5 * [b + sgn(b) * sqrt(b*b - 4ac)].
This algorithm is particularly good at yielding accurate results when
a and/or c are small values.
*/
BOOL FindQuadRoots (
float a, /* input - coefficient a of quadratic equation */
float b, /* input - coefficient b of quadratic equation */
float c, /* input - coefficient c of quadratic equation */
float *x1, /* output - first root computed */
float *x2) /* output - second root computed */
{
float d; /* root algorithm variable */
BOOL root_stat; /* status of root computations */
d = b * b - 4 * a *c;
if (d < 0)
return FAILURE;
else
{
d = (float)sqrt((double)d);
/* make the result of sqrt the sign of b */
if (b < 0 )
d = -d;
d = -0.5 * (b + d);
/* solve for the roots of the quadratic */
/* if both root computations will yield divide by zero ... forget it! */
if ( (a == 0) && (d == 0) )
return FAILURE;
root_stat = SUCCESS;
/* compute first root if denominator a is not zero */
if (a == 0)
root_stat = FAILURE1;
else
*x1 = d / a;
/* compute second root if denominator d is not zero */
if (d == 0)
root_stat = FAILURE2;
else
*x2 = c / d;
return root_stat;
}
}
/* Given an abscissa (x) location, computes the corresponding cubic spline
ordinate (y) value.
*/
void ComputeY (
unsigned int i, /* input - array index */
float *x, /* input - array of x values */
float x_value, /* input - x value at which to solve for y */
float *y, /* input - array of y values */
float *sec_deriv, /* input - array of second derivatives of each data interval */
float *y_value) /* output - address of y extreme value at x */
{
float A, B, C, D; /* cubic spline coefficients */
float ftemp; /* temporary float */
/* compute the standard cubic spline coefficients */
A = (x[i + 1] - x_value) / (x[i + 1] - x[i]);
B = 1 - A;
ftemp = (float) pow((double)(x[i + 1] - x[i]), 2.0) / 6.0;
C = (A * A * A - A) * ftemp;
D = (B * B * B - B) * ftemp;
/* compute the ordinate value at the abscissa location */
*y_value = A * y[i] + B * y[i + 1] + C * sec_deriv[i] + D *
sec_deriv[i + 1];
return;
}
