Mike J. Courtney

When collecting or computing discrete data, particularly time-varying data, the resulting resolution is often limited by practical influences such as computational burden, A/D sampling rates, storage limitations, and transmit rates. This can present a challenge when the minimum and/or maximum values within the continuous function describing the data is desired. 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.

Determination of these minima and maxima, which are collectively termed extrema, requires some form of interpolation between the known discrete data points. However, if the interpolation resolution is too course, the true extrema might be overlooked and if the resolution is too fine, then the computations would take entirely too long . There are a number of algorithms for more efficiently finding the extrema of a data set [1]. These algorithms typically consist of bracketed searches and are often iterative.

An alternative approach is presented here that I have derived from the well-established equation for cubic spline interpolation. This approach performs a direct, non-iterative determination of the extrema of a discrete data set.

Presented is an overview of the algorithm approach and derivation, the software implementation, and test data results. These data illustrate the performance of the program and validate algorithm accuracy via comparison to expected results, including examples from the worlds of physics and electronics control theory. The software implementation consists of an input/output program and the algorithmic routines.

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, if desired.

The method uses all the given data to compute second derivatives at each point, which are also called knots. Then the space between each knot, termed 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, accurate, and its implementation is quite robust.

The goal of cubic spline interpolation is to derive an interpolation formula such that 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 yi = y(xi), the linear interpolation in the interval between xj and xj+1 can be expressed as

where " denotes the second derivative and

A = (xj+1 - x ) / (xj+1 - xj), B = 1 - A, C = ( (A3 - A)(xj+1 - xj)2)/6, and D = ((B3 - B)(xj+1 - xj)2)/6.

The first derivative of (1) is denoted as dy/dx and is solved as

dy/dx = (yi+1 - yi)/(xi+1 - xi) - [(3A2 - 1)(xi+1 - xi)y"i]/6 + [(3B2 - 1)(xi+1 - xi)y"i+1]/6. (2)

The derivation of these equations is available in a number of text [2]. However, if we take this further and set (2) to zero, then a new equation can be derived that represents the maximization and minimization of (1). This new equation will allow identification of the points at which y remains constant with respect to finite changes in x.

Finally, equation (2) can be expressed in the quadratic form ax2 + bx + c = 0, such that x can be solved for to yield the cubic spline extrema equation where

a = 3(y"i+1 - y"i), (3)

b = 6(xi+1 y"i - xiy"i+1), (4)

c = 6yi+1 - 6yi - (2x2i+1 - x2i + 2xixi+1)y"i - (x2i+1 - 2x2i - 2xixi+1)y"i+1. (5)

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.

The data input and output program is called cubic_io.c and is shown in Listing 1. Its purpose is to read and parse an input ascii data file, call FindCubicExtrema(), and then print 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. An important criterion placed upon the input data is that each x value must be greater than the previous one. Because it is such a trivial task to ensure this prior to calling FindCubicExtrema(), it is a task left to the user. However, if subjected to such erroneous data, FindCubicExtrema() has ways to deal with it that are discussed later.

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. Every allocation is asserted to ensure that the operations were successful. If an assertion “hits”, then the program exits, and the failed array name and program line number are automatically printed. If the allocations are successful, 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 (see Listing 2).

Because we are not using a doubly linked list that could provide pointers to the previous data, we will need some other way to get back to the first structure in the list. This is accomplished by saving its address in the pointer variable first_extr prior to calling FindCubicExtrema().

FindCubicExtrema() is then called and passed pointers to the number of input data points, the input x and y arrays, and the first output extremum structure. If no extremum were found, then the function returns a FAILURE and the user is notified. If the function returned 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, as opposed to just the absolute extrema, because no assumptions are made about the user’s application. If the intent were to compute the period of the input data, then all extrema might be desired to compute an average period. If the rate of decay of an underdamped oscillatory function were to be examined, all maxima would be needed. 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.

The primary routine which computes the data extrema is FindCubicExtrema() and is presented in Listing 3. An overview of the processing that takes place in this routine is as follows:

1. ComputeSecDerivs() is called to compute the second derivatives of each interval;

2. The cubic extrema quadratic coefficients are calculated and FindQuadRoots() is called to compute the roots of the quadratic equation;

3. It is determined if the roots, which are candidate extremum abscissa, lie within the current interval;

4. For each valid extremum abscissa, ComputeY() is called to determine the corresponding ordinate value.

If we examine the derivative form of the cubic spline equation (2), we can readily observe that we have 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 thru aN1x1 + aN2x2 + ... + aNNxN = bN.

These equations may be represented in matrix form as Ax = b, where A represents the matrix of a’s. 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 b’s are represented as right[] to indicate the right hand side of the equations. The second derivative variables x are denoted as sec_deriv[]. The latter is what we are solving for in this routine.

Observing the ComputeSecDerivs() code implementation, the computation of the second derivatives involves division by the diagonals, which are computed from the abscissa data x[]. Earlier we stated that a pre-existing condition of the data was that the abscissa values be increasing. This also implies that the data be unique. It intuitively makes sense that this should be the case for the problem we have posed. If the abscissa did not change but the ordinate did, then the data would represent an impulse function and the corresponding derivative will be infinity.

If these abscissa data were non-unique, 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. An assertion approach is used instead of an “if ... then return FAILURE” approach because this is non-recoverable error and, according to the rules we have established, this condition should theoretically never occur. However, if it does, then something has gone awry earlier in allowing this data to get through earlier defenses.

Because the second derivatives at the knots are now available, the cubic extrema quadratic coefficients are easily solved using the equations (3) thru (5). 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 [3] that avoids associated accuracy errors.

Note that in the software implementation of FindQuadRoots() there are 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 an 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 root 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. Note that assertions are not used because these conditions do not indicate a non-recoverable and unexpected error situation that is not recoverable.

If a status other than FAILURE has been returned from FindQuadRoots(), then at least one of the roots is a valid abscissa candidate. Each is then tested regarding their return status and if accepted for further testing, bounds checking is performed. This simply consists of determining if the root lies within the current interval. If it does, then it is a valid root and an abscissa location at which an extremum lies!

Since we have a valid extremum abscissa, we want to save it in the linked list. To see if we are at the first root in the list, the flag first_flag is tested. An alternative might be to check the next pointer, which has not yet been set to NULL when at the first structure. However, this pre-supposes that the next pointer is non-NULL upon creation. This is a dangerous assumption that may be influenced by platform and even a debugger. Therefore, the flag is instead tested against a known state over which we have control.

If at the beginning of the list, which does not yet have any links, there is nothing unique that must occur except for the setting of first_flag to FALSE for subsequent passes. If we are beyond the first extremum, then we need to add onto the chain by allocating another structure linked from the next pointer and must also set the current structure to the next structure. In either case, the end of the chain is then marked by setting the structure’s next pointer to NULL and the extremum abscissa value is stored in the structure element x.

A perhaps unapparent reason for returning the specific FAILURE1 or FAILURE2 status from FindQuadRoots() is to avoid 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.

If a valid abscissa has been found, it is then utilized to compute the corresponding ordinate value with ComputeY(). Also required is the input data ordinates and the second derivatives, which have already been computed. Thus all the hard work has already been performed and this computation is easy using the standard equation for a cubic spline (1).

If one or more roots were found in FindCubicExtrema(), SUCCESS is returned to main() so that no interrogation of the list is required unless valid data exists. Furthermore, since the first structure was allocated in main(), there would be no other way to determine if the structure contained valid data, outside of passing back the number of extrema found. 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.

It is improper to label software as complete unless it fulfills its purpose with accurate and consistent behavior. Ensuring this behavior requires that the software be validated utilizing input data that yields calculable and thus comparable outputs. These different data sets represent test cases that provide performance verification and “stress” testing. Many such test cases were developed to aid in the testing, debugging, and validation of the cubic spline extrema algorithm and its software implementation. Presented here are some of these input data, the expected and actual software outputs, and their corresponding errors.

The first data set presented in a simple one consisting of the 4 data points plotted below, with a polynomial interpolating curve overlaid for visual aid. Due to data symmetry, it is obvious that a maximum is expected at an x value of 0.0 and a y value of approximately 1.1. This data set, aside from resembling St. Louis’s Gateway Arch, permits the validation of several important attributes. These are that the program can handle negative input data and that it can solve for an extremum at an abscissa of zero. The program was run utilizing these data and the extremum computed were as expected.

Expected: Actual:

x = 0.0, y~1.1 x = 0.00, y = 1.15

The last 3 of the 4 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 2 data points can not yield an extremum — only a straight line. At least 3 are required, therefore representing a minimal data set. This “stress” test also illustrates that nothing in the algorithm or software precludes an extremum in the first or last interval from being located. 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%, which is quite satisfactory given the small amount of data provided to the algorithm.

The next data set presented is an extension of the first in that it yields a single maximum. However, this 7 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 the figure was hand-calculated based upon a projectile shot at 1500 ft/sec at an initial angle of +45 degrees and a position sampling at 10 second intervals until ground impact.

h

The standard trajectory equations for the maximum height and corresponding distance were then used to hand-calculate the expected extremum.

Expected: Actual:

x = 34,937.89, y = 17,468.94 x = 34,896.04, y = 17,469.07

The x error was only about 1/10% and the y error was negligible.

Cubic splines operate best when provided exact data. But what if the data were only 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 20th feet. For example, 21213.2 feet becomes 21200.0 feet. This provided a good simulation of the errors that might be expected if the GPS (Global Positioning System) 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 1/10%.

The next data set is one 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 shown in the following figure.

As illustrated, this data set should yield both maxima and minima. An interpolating curve is not shown in this figure because even a single 5th 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 calculated. Five extrema are found using this function-dependent calculation. By supplying the data set to the cubic extrema algorithm, the following results were yielded compared to the expected:

Expected: Actual:

t = 0.135, y(t) = -1.986 t = none, y(t) = none

t = 0.920, y(t) = 1.547 t = 0.915, y(t) = 1.551

t = 1.706, y(t) = 0.751 t = 1.707, y(t) = 0.751

t = 2.491, y(t) = 1.114 t = 2.492, y(t) = 1.114

t = 3.277, y(t) = 0.949 t = 3.279, y(t) = 0.949

The first extremum was not found, although the cause is obvious from observation of the figure. 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 that were detected, the largest extremum error noted in the time axis was only 1/2 % and the maximum y(t) error was about 1/4 %.

The last data presented provides 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 8 point data set depicted in the following figure are samples of the function f(x) = x3 - 2x2 + x + 1.

Note that the 4th interval is quite large compared to the others. Taking the derivative of this function shows that there is both a maximum and a minimum in this interval. If you think about it, there is no way to have 3 extrema in an interval and thus 2 is the maximum. The combination of irregular sampling and two subtle extrema make this an extremely challenging and uncommon situation. The function’s expected and resulting data are:

Expected: Actual:

x = 0.333, f(x) = 1.148 x = 0.318, f(x) = 1.147

x = 1.000, f(x) = 1.000 x = 1.030, f(x) = 1.007

These data indicate abscissa errors of less than 5% and ordinate errors of less than 1%. 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 is therefore invalidated by the bounds checking. Therefore, in such situations, the more data the better.

Several additional test cases were utilized 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. Also tested was a data set that has decreasing abscissa values, which is a violation of the pre-existing condition. Using this data showed that while no extrema were calculable with this data, the software did not fault either. Additionally, a data set was utilized which itself contains an extremum. This was important to ensure that the root computations and bounds checking logic would not preclude an extremum which coincides with an interval knot.

The cubic spline extrema algorithm has been shown to be an effective method of determining the relative extrema of a given data set. Its software implementation has been shown to be quite thorough and robust. Several data sets have been presented that illustrate a number of important attributes. These can be summarized as follows:

[1] Press, Williams. Numerical Recipes in C, 2nd ed., Cambridge University Press, pp 394 - 455.

[2] ibid, pp 113-114.

[3] ibid, pp 183 - 184.