Spline interpolation

Markus NentwigMay 11, 20142 comments

A cookbook recipe for segmented y=f(x) 3rd-order polynomial interpolation based on arbitrary input data.
Includes Octave/Matlab design script and Verilog implementation example.


Keywords: Spline, interpolation, function modeling, fixed point approximation, data fitting, Matlab, RTL, Verilog

Introduction

Splines describe a smooth function with a small number of parameters.
They are well-known for example from vector drawing programs, or to define a "natural" movement path through given points in computer animation.
They are an engineer's tool, apparently originating from the shipbuilding industry, with widespread use documented from world war II to make "backup copies" of aircraft design templates that did not fit into a bomb shelter.

Splines are convenient to model arbitrary continuous functions, as fixed-point implementation is efficient and straightforward. In DSP applications, a spline's continuous first derivative helps to confine interpolation error close to frequencies of the signal, which is often desirable (i.e. to limit audio aliasing to higher frequencies where the ear is more sensitive). For example, Catmull-Rom splines are commonly used for low-cost audio resampling.

This article provides a recipe to construct spline coefficients from arbitrary scattered points via a least-squares solver.
For practical design problems this allows a fluid "cooking-style" design process, where the result is shaped by adding (or removing) input data to control the approximation error.


Figure 1: Spline approximation of a function, possibly from noisy data


The splines discussed in this article consist of multiple 3rd-order polynomial sections.
Transitions between segments are continuous both in value and in the first derivative, resulting in a function that is for most practical purposes fairly "smooth".

Note that the definition of "spline" can be quite general.
Here, continuity in the second derivative is not enforced. This results in a more accurate function approximation, at the expense of "smoothness".

Examples

The example in figure 2 is a logarithmic volume control curve from a DIY FPGA project. Using four segments, the spline approximation is virtually indistinguishable from the original exponential function. The approximation error is well below 1 percent (-40 dB).


Figure 2: A four-segment spline modeling a 40 dB-range volume curve


A second example (figure 3) shows a spline fit to noisy and non-gridded data.
It originates from the same project, where the target waveform represents three cycles of a Hammond organ drawbar signal, with noise added both to sampled values (y axis) and sampling time instants (x axis).
As can be seen, the design process is robust towards irregularities in the data.


Figure 3: 16-segment spline through a very noisy Hammond organd drawbar sample

The math

For segmentation, the input variable is split into an integer part and a fractional part (figure 4).


Figure 2: A four-segment spline modeling a 40 dB-range volume curve


For a number of segments that is a power of two, the separation is done by taking the higher- and lower bits.
The Verilog implementation below shows the details for signed numbers.
Each of the segments is described by a polynomial y(x) = a0+a1x+a2x2+a3x3 with x in a domain interval -1..1 covering the segment (see "polynomial variable" in figure 4).

Polynomials are defined symmetrically over a [-1..1] range to improve numerical accuracy (because xn=1 at x=1 for any n).
Note that a segment drawn as k..k+1 appears in the polynomial evaluation with twice the range [-1..1].


Within a segment, four basis functions p1.. p4 are defined that have zero value and zero slope at x = +/-1, with the exception of

  • p1 has value 1 at x = -1
  • p2 has value 1 at x = 1
  • p3 has slope 1 at x = -1
  • p4 has slope 1 at x = -1

Figure 5: Basis functions


For each support point, a first scaling factor is assigned to a pair of basis functions p1, p2 on the left and right side of the support point, respectively: The scaling factor controls the value of the spline at the support point, but does not affect the slope (because the slope of both p1, p2 is zero at the support point). This is shown in the upper half of Figure 6.
In a similar manner, a second scaling factor is assigned to p3, p4 to set the slope, without affecting the value at the support point (lower half of Figure 6).


Figure 6: Scaled basis functions enforcing continuity through a support point


The scaling factors are then found by minimizing the difference between input data and approximation in a least-squares problem.
Finally, the polynomial coefficients for each spline segment are determined by summing up all four scaled basis function polynomials.
The derivation for the basis functions in Maxima CAS can be found here.

A generalization of this idea is known as B-Splines. The basis functions p1 to p4 could be combined into an equivalent B-spline that spans three segments.
The discussed approach is slightly different, as it uses the slope as an explicit parameter.
Number of support points

To span n segments, n+1 support points are needed.
Relevant for all practical purposes is the number of segments, not support points, as each segment defines one set of polynomial coefficients.
For example, a fixed-point implementation may use four MSBs to select one of 16 segments. The equation system in the design process will use 17 support points.
In the provided Matlab program, the outermost support points do not need any special consideration with regard to continuity: In an interval without input data, any basis function becomes invisible to the least-squares cost function.

Octave (Matlab) script

All the little details are taken care of in the Octave / Matlab program, which can be downloaded here.
Several examples to generate data are included. A spline is generated, evaluated and plotted.
Finally, fixed-point coefficients are written to a file that is imported by the Verilog RTL implementation.

RTL implementation

This Verilog implementation demonstrates the evaluation of the spline in RTL. It was written with the 18-bit multipliers on Spartan-6 FPGAs in mind.

Input variable

The signed input is first converted to offset binary by inverting the sign bit.
The four MSBs become the index into the coefficient table. The LSBs are then converted back in the same manner from offset binary to two's complement signed.


Figure 7: Fixed-point handling of input variable
Multiply-accumulate

With four MSBs stripped off, the LSB part of the input variable contains 13 fractional bits. After multiplication with the LSB part, arithmetic right shift by 13 bits restores the decimal point. Prior to that, 0.5 LSB is added for correct midpoint rounding.

Note that the removed MSBs provide numeric headroom for the MAC-operation.
In some cases (i.e. only two segments) it may be insufficient and cause internal overflow in the MAC.
If this happens, a straightforward solution is to divide all coefficients by 2 at design time, then simply double the result.

Figure 8: Horner scheme evaluation


In the Verilog code, the four operations of the Horner scheme evaluation are written out explicitly for clarity.
Alternatively, a MAC could be used, with the feedback path opened during the first cycle to load c3 (i.e. "CLR" input).

Simulation result

The simulation snapshot in figure 9 covers the whole 18-bit input range of input values.
The segment structure is clearly visible in MSB, LSB and the selected coefficients (see figure 7), and the output signal looks as expected.


Figure 9: Simulation result for "Hammond drawbar" spline

Pipelined RTL implementation

A second RTL version ( download) is a pipelined implementation of the four-segment volume curve example from Figure 2. The evaluation of four independent input variable values is interleaved.
The pipelining generally allows higher clock frequencies by introducing additional delay.

Numerical accuracy

The implementations are not bit-exact because of intermediate rounding after every multiplication. A typical error is one LSB in the output.

Conclusion

This article presents a design method to construct an interpolation function in fixed-point RTL from arbitrary data.
A spline is derived via a least-squares fit. The spline is then evaluated using the Horner scheme.
This work was originally done for an open-source FPGA synth project. If it gets reused (amplifier model, anybody?), I'd be happy to hear about it.

Download

The whole project archive, including the links from the text, can be downloaded here.


Previous post by Markus Nentwig:
   Signed serial-/parallel multiplication
Next post by Markus Nentwig:
   A poor man's Simulink


Comments:

[ - ]
Comment by maksim4August 16, 2014
You calculate the coefficients by Octave / Matlab. Did considerate calculation of the coefficients for given input data by the FPGA itself?
[ - ]
Comment by mnentwigAugust 17, 2014
Hi,

for a least-squares approach as in the Matlab script, it usually boils down to a multiply-accumulate operation between the input data and a (precalculated) row of the matrix inverse. There are some variants (i.e. Lagrange or Catmull-Rom) where the coefficients are polynomials of the support points. It depends on the problem.
Generally, I'd try to avoid doing anything in RTL if it's not performance critical. It sounds more like a job for a soft-core CPU.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.

Sign up
or Sign in