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 3

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”.

Here, continuity in the

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)

As can be seen, the design process is robust towards irregularities in the data.

Figure 3: 16-segment spline through a

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

Polynomials are defined symmetrically over a [-1..1] range to improve numerical accuracy (because *x*^{n}=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].

Note that a segment drawn as

Within a segment, four basis functions *p _{1}.. p_{4}* are defined that have zero value and zero slope at

*p*has value 1 at_{1}*x = -1**p*has value 1 at_{2}*x = 1**p*has slope 1 at_{3}*x = -1**p*has slope 1 at_{4}*x = -1*

Figure 5: Basis functions

For each support point, a first scaling factor is assigned to a pair of basis functions

In a similar manner, a second scaling factor is assigned to

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.

A generalization of this idea is known as B-Splines. The basis functions *p*_{1} to *p*_{4} 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.

The discussed approach is slightly different, as it uses the slope as an explicit parameter.

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.

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.

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.

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

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.

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 c

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

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.

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

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.