Cubic Spline Interpolation

Functions
RISCV_DSP_ATTRIBUTE void	riscv_spline_f32 (riscv_spline_instance_f32 *S, const float32_t *xq, float32_t *pDst, uint32_t blockSize)
	Processing function for the floating-point cubic spline interpolation. More...
RISCV_DSP_ATTRIBUTE void	riscv_spline_init_f32 (riscv_spline_instance_f32 *S, riscv_spline_type type, const float32_t *x, const float32_t *y, uint32_t n, float32_t *coeffs, float32_t *tempBuffer)
	Initialization function for the floating-point cubic spline interpolation. More...

Detailed Description

Spline interpolation is a method of interpolation where the interpolant is a piecewise-defined polynomial called "spline".

Introduction

Given a function f defined on the interval [a,b], a set of n nodes x(i) where a=x(1)<x(2)<...<x(n)=b and a set of n values y(i) = f(x(i)), a cubic spline interpolant S(x) is defined as:

        S1(x)       x(1) < x < x(2)
S(x) =   ...         
        Sn-1(x)   x(n-1) < x < x(n)

where

Si(x) = a_i+b_i(x-xi)+c_i(x-xi)^2+d_i(x-xi)^3    i=1, ..., n-1

Algorithm

Having defined h(i) = x(i+1) - x(i)

h(i-1)c(i-1)+2[h(i-1)+h(i)]c(i)+h(i)c(i+1) = 3/h(i)*[a(i+1)-a(i)]-3/h(i-1)*[a(i)-a(i-1)]    i=2, ..., n-1

It is possible to write the previous conditions in matrix form (Ax=B). In order to solve the system two boundary conidtions are needed.

• Natural spline: S1''(x1)=2*c(1)=0 ; Sn''(xn)=2*c(n)=0 In matrix form:

|  1        0         0  ...    0         0           0     ||  c(1)  | |                        0                        |
| h(0) 2[h(0)+h(1)] h(1) ...    0         0           0     ||  c(2)  | |      3/h(2)*[a(3)-a(2)]-3/h(1)*[a(2)-a(1)]      |
| ...      ...       ... ...   ...       ...         ...    ||  ...   |=|                       ...                       |
|  0        0         0  ... h(n-2) 2[h(n-2)+h(n-1)] h(n-1) || c(n-1) | | 3/h(n-1)*[a(n)-a(n-1)]-3/h(n-2)*[a(n-1)-a(n-2)] |
|  0        0         0  ...    0         0           1     ||  c(n)  | |                        0                        |

• Parabolic runout spline: S1''(x1)=2*c(1)=S2''(x2)=2*c(2) ; Sn-1''(xn-1)=2*c(n-1)=Sn''(xn)=2*c(n) In matrix form:

|  1       -1         0  ...    0         0           0     ||  c(1)  | |                        0                        |
| h(0) 2[h(0)+h(1)] h(1) ...    0         0           0     ||  c(2)  | |      3/h(2)*[a(3)-a(2)]-3/h(1)*[a(2)-a(1)]      |
| ...      ...       ... ...   ...       ...         ...    ||  ...   |=|                       ...                       |
|  0        0         0  ... h(n-2) 2[h(n-2)+h(n-1)] h(n-1) || c(n-1) | | 3/h(n-1)*[a(n)-a(n-1)]-3/h(n-2)*[a(n-1)-a(n-2)] |
|  0        0         0  ...    0        -1           1     ||  c(n)  | |                        0                        |

A is a tridiagonal matrix (a band matrix of bandwidth 3) of size N=n+1. The factorization algorithms (A=LU) can be simplified considerably because a large number of zeros appear in regular patterns. The Crout method has been used: 1) Solve LZ=B

u(1,2) = A(1,2)/A(1,1)
z(1)   = B(1)/l(11)
 
FOR i=2, ..., N-1
  l(i,i)   = A(i,i)-A(i,i-1)u(i-1,i)
  u(i,i+1) = a(i,i+1)/l(i,i)
  z(i)     = [B(i)-A(i,i-1)z(i-1)]/l(i,i)
 
l(N,N) = A(N,N)-A(N,N-1)u(N-1,N)
z(N)   = [B(N)-A(N,N-1)z(N-1)]/l(N,N)

2) Solve UX=Z

c(N)=z(N)
 
FOR i=N-1, ..., 1
  c(i)=z(i)-u(i,i+1)c(i+1) 

c(i) for i=1, ..., n-1 are needed to compute the n-1 polynomials. b(i) and d(i) are computed as:

• b(i) = [y(i+1)-y(i)]/h(i)-h(i)*[c(i+1)+2*c(i)]/3
• d(i) = [c(i+1)-c(i)]/[3*h(i)] Moreover, a(i)=y(i).

Behaviour outside the given intervals

It is possible to compute the interpolated vector for x values outside the input range (xq<x(1); xq>x(n)). The coefficients used to compute the y values for xq<x(1) are going to be the ones used for the first interval, while for xq>x(n) the coefficients used for the last interval.

Initialization function

The initialization function takes as input two arrays that the user has to allocate: coeffs will contain the b, c, and d coefficients for the (n-1) intervals (n is the number of known points), hence its size must be 3*(n-1); tempBuffer is temporally used for internal computations and its size is n+n-1.

The x input array must be strictly sorted in ascending order and it must not contain twice the same value (x(i)<x(i+1)).

Function Documentation

riscv_spline_f32()

RISCV_DSP_ATTRIBUTE void riscv_spline_f32	(	riscv_spline_instance_f32 *	S,
		const float32_t *	xq,
		float32_t *	pDst,
		uint32_t	blockSize
	)

Processing function for the floating-point cubic spline interpolation.

Parameters

[in]	S	points to an instance of the floating-point spline structure.
[in]	xq	points to the x values of the interpolated data points.
[out]	pDst	points to the block of output data.
[in]	blockSize	number of samples of output data.

riscv_spline_init_f32()

RISCV_DSP_ATTRIBUTE void riscv_spline_init_f32	(	riscv_spline_instance_f32 *	S,
		riscv_spline_type	type,
		const float32_t *	x,
		const float32_t *	y,
		uint32_t	n,
		float32_t *	coeffs,
		float32_t *	tempBuffer
	)

Initialization function for the floating-point cubic spline interpolation.

Parameters

[in,out]	S	points to an instance of the floating-point spline structure.
[in]	type	type of cubic spline interpolation (boundary conditions)
[in]	x	points to the x values of the known data points.
[in]	y	points to the y values of the known data points.
[in]	n	number of known data points.
[in]	coeffs	coefficients array for b, c, and d
[in]	tempBuffer	buffer array for internal computations