1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_spline_interp_init_f32.c
4 * Description: Floating-point cubic spline initialization function
5 *
6 * $Date: 23 April 2021
7 * $Revision: V1.9.0
8 *
9 * Target Processor: Cortex-M and Cortex-A cores
10 * -------------------------------------------------------------------- */
11 /*
12 * Copyright (C) 2010-2021 ARM Limited or its affiliates. All rights reserved.
13 *
14 * SPDX-License-Identifier: Apache-2.0
15 *
16 * Licensed under the Apache License, Version 2.0 (the License); you may
17 * not use this file except in compliance with the License.
18 * You may obtain a copy of the License at
19 *
20 * www.apache.org/licenses/LICENSE-2.0
21 *
22 * Unless required by applicable law or agreed to in writing, software
23 * distributed under the License is distributed on an AS IS BASIS, WITHOUT
24 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
25 * See the License for the specific language governing permissions and
26 * limitations under the License.
27 */
28
29 #include "dsp/interpolation_functions.h"
30
31 /**
32 @ingroup groupInterpolation
33 */
34
35 /**
36 @addtogroup SplineInterpolate
37 @{
38
39 @par Initialization function
40
41 The initialization function takes as input two arrays that the user has to allocate:
42 <code>coeffs</code> will contain the b, c, and d coefficients for the (n-1) intervals
43 (n is the number of known points), hence its size must be 3*(n-1); <code>tempBuffer</code>
44 is temporally used for internal computations and its size is n+n-1.
45
46 @par
47
48 The x input array must be strictly sorted in ascending order and it must
49 not contain twice the same value (x(i)<x(i+1)).
50
51 */
52
53 /**
54 * @brief Initialization function for the floating-point cubic spline interpolation.
55 * @param[in,out] S points to an instance of the floating-point spline structure.
56 * @param[in] type type of cubic spline interpolation (boundary conditions)
57 * @param[in] x points to the x values of the known data points.
58 * @param[in] y points to the y values of the known data points.
59 * @param[in] n number of known data points.
60 * @param[in] coeffs coefficients array for b, c, and d
61 * @param[in] tempBuffer buffer array for internal computations
62 *
63 */
64
arm_spline_init_f32(arm_spline_instance_f32 * S,arm_spline_type type,const float32_t * x,const float32_t * y,uint32_t n,float32_t * coeffs,float32_t * tempBuffer)65 void arm_spline_init_f32(
66 arm_spline_instance_f32 * S,
67 arm_spline_type type,
68 const float32_t * x,
69 const float32_t * y,
70 uint32_t n,
71 float32_t * coeffs,
72 float32_t * tempBuffer)
73 {
74 /*** COEFFICIENTS COMPUTATION ***/
75 /* Type (boundary conditions):
76 - Natural spline ( S1''(x1) = 0 ; Sn''(xn) = 0 )
77 - Parabolic runout spline ( S1''(x1) = S2''(x2) ; Sn-1''(xn-1) = Sn''(xn) ) */
78
79 /* (n-1)-long buffers for b, c, and d coefficients */
80 float32_t * b = coeffs;
81 float32_t * c = coeffs+(n-1);
82 float32_t * d = coeffs+(2*(n-1));
83
84 float32_t * u = tempBuffer; /* (n-1)-long scratch buffer for u elements */
85 float32_t * z = tempBuffer+(n-1); /* n-long scratch buffer for z elements */
86
87 float32_t hi, hm1; /* h(i) and h(i-1) */
88 float32_t Bi; /* B(i), i-th element of matrix B=LZ */
89 float32_t li; /* l(i), i-th element of matrix L */
90 float32_t cp1; /* Temporary value for c(i+1) */
91
92 int32_t i; /* Loop counter */
93
94 S->x = x;
95 S->y = y;
96 S->n_x = n;
97
98 /* == Solve LZ=B to obtain z(i) and u(i) == */
99
100 /* -- Row 1 -- */
101 /* B(0) = 0, not computed */
102 /* u(1,2) = a(1,2)/a(1,1) = a(1,2) */
103 if(type == ARM_SPLINE_NATURAL)
104 u[0] = 0; /* a(1,2) = 0 */
105 else if(type == ARM_SPLINE_PARABOLIC_RUNOUT)
106 u[0] = -1; /* a(1,2) = -1 */
107
108 z[0] = 0; /* z(1) = B(1)/a(1,1) = 0 always */
109
110 /* -- Rows 2 to N-1 (N=n+1) -- */
111 hm1 = x[1] - x[0]; /* Initialize h(i-1) = h(1) = x(2)-x(1) */
112
113 for (i=1; i<(int32_t)n-1; i++)
114 {
115 /* Compute B(i) */
116 hi = x[i+1]-x[i];
117 Bi = 3*(y[i+1]-y[i])/hi - 3*(y[i]-y[i-1])/hm1;
118
119 /* l(i) = a(i)-a(i,i-1)*u(i-1) = 2[h(i-1)+h(i)]-h(i-1)*u(i-1) */
120 li = 2*(hi+hm1) - hm1*u[i-1];
121
122 /* u(i) = a(i,i+1)/l(i) = h(i)/l(i) */
123 u[i] = hi/li;
124
125 /* z(i) = [B(i)-h(i-1)*z(i-1)]/l(i) */
126 z[i] = (Bi-hm1*z[i-1])/li;
127
128 /* Update h(i-1) for next iteration */
129 hm1 = hi;
130 }
131
132 /* -- Row N -- */
133 /* l(N) = a(N,N)-a(N,N-1)u(N-1) */
134 /* z(N) = [-a(N,N-1)z(N-1)]/l(N) */
135 if(type == ARM_SPLINE_NATURAL)
136 {
137 /* li = 1; a(N,N) = 1; a(N,N-1) = 0 */
138 z[n-1] = 0; /* a(N,N-1) = 0 */
139 }
140 else if(type == ARM_SPLINE_PARABOLIC_RUNOUT)
141 {
142 li = 1+u[n-2]; /* a(N,N) = 1; a(N,N-1) = -1 */
143 z[n-1] = z[n-2]/li; /* a(N,N-1) = -1 */
144 }
145
146 /* == Solve UX = Z to obtain c(i) and */
147 /* compute b(i) and d(i) from c(i) == */
148
149 cp1 = z[n-1]; /* Initialize c(i+1) = c(N) = z(N) */
150
151 for (i=n-2; i>=0; i--)
152 {
153 /* c(i) = z(i)-u(i+1)c(i+1) */
154 c[i] = z[i]-u[i]*cp1;
155
156 hi = x[i+1]-x[i];
157 /* b(i) = [y(i+1)-y(i)]/h(i)-h(i)*[c(i+1)+2*c(i)]/3 */
158 b[i] = (y[i+1]-y[i])/hi-hi*(cp1+2*c[i])/3;
159
160 /* d(i) = [c(i+1)-c(i)]/[3*h(i)] */
161 d[i] = (cp1-c[i])/(3*hi);
162
163 /* Update c(i+1) for next iteration */
164 cp1 = c[i];
165 }
166
167 /* == Finally, store the coefficients in the instance == */
168
169 S->coeffs = coeffs;
170 }
171
172 /**
173 @} end of SplineInterpolate group
174 */
175
176