はじめに
補間法の一つスプライン補間をC言語で書いてみました。
「何をいまさら」ですが、まあ見てください。
奥村晴彦先生の「C言語による標準アルゴリズム事典」を参考にしました。
方針
spline関数は表を読み、データーの繋目の条件をあらかじめ計算して置き、あとから呼び出されたときにその値を使って補間値を求めます。奥村先生の本では、一つのファイルに表を作る関数と、計算する関数を同居させました。
この本はアルゴリズムを説明するためにこうしてますが、IoT等組み込み分野では、この実装は不利です。
- 解析した結果は、初めのデーターとともにROMに配置したい
- 表の作成は、組み込み機ではなく、開発機で行いたい
- 複数の表を、交互に使いたい
これらを実現できるように書いてみます。
表（曲率データー）の作成と、実機での計算を分離したのが最大の特徴でが、同じマシン上でももちろん結構です。
表を作成するプログラム
makespline.c
/*
* makespline.c:
* Make spline table
* Coded by fanfan (S.Makino)
* Copyright (c) 2021 Agri-hitech LLC.
*/
typedef struct {
double in;
double out;
double diff;
} TBL;
// ここに表を書く
TBL data[] = {
{0,0,0},
{1,1,0},
{2,2,0},
{3,4,0},
{4,8,0},
{5,8,0},
{6,4,0},
{7,2,0},
{8,1,0},
{9,0,0},
};
/***********************************************************
maketable.c -- スプライン補間データー作成
***********************************************************/
void maketable(TBL *T, int size)
{
int i;
double t;
double h[size], d[size];
T[0].diff = 0; T[size - 1].diff = 0; /* 両端点での y''(x) / 6 */
for (i = 0; i < size - 1; i++) {
h[i ] = T[i + 1].in - T[i].in;
d[i + 1] = (T[i + 1].out - T[i].out) / h[i];
}
T[1].diff = d[2] - d[1] - h[0] * T[0].diff;
d[1] = 2 * (T[2].in - T[0].in);
for (i = 1; i < size - 2; i++) {
t = h[i] / d[i];
T[i + 1].diff = d[i + 2] - d[i + 1] - T[i].diff * t;
d[i + 1] = 2 * (T[i + 2].in - T[i].in) - h[i] * t;
}
T[size - 2].diff -= h[size - 2] * T[size - 1].diff;
for (i = size - 2; i > 0; i--)
T[i].diff = (T[i].diff - h[i] * T[i + 1].diff) / d[i];
}
#define Nd (sizeof(data)/sizeof(TBL))
#include <stdio.h>
int main(){
int i;
maketable((void *)data,Nd);
printf("// Auto generated by mkspline.c\n");
printf("// DON'T Edit !!\n");
printf("const TBL data[] = {\n");
for (i = 0; i < Nd;i++) {
printf("\t{%.16g,%.16g,%.16g},\n",
data[i].in,data[i].out,data[i].diff);
}
printf("};\n");
return 0;
}
// ここに表を書く
TBL data[] = {
以下にに、データー点を書いてください。
{入力, 出力, 0},
と列挙します。
等間隔でなくても、入力が昇順になっていれば結構です。
３番目はこのプログラムが踏みつぶしますので、普通0を書いておきます。
makesplineを実行する
gcc makespline.c -o makespline
./makespline > data.c
data.cを見てみましょう。
data.c
// Auto generated by makespline.c
// DON'T Edit !!
const TBL data[] = {
{0,0,0},
{1,1,-0.0188679},
{2,2,0.0754717},
{3,4,0.716981},
{4,8,-0.943396},
{5,8,-0.943396},
{6,4,0.716981},
{7,2,0.0754717},
{8,1,-0.0188679},
{9,0,0},
};
プログラムに書いた表が変換されてます。
constがついていますので、ROMに配置され、3番目のデーターに計算結果が入ってます。
スプライン本体
スプライン関数本体はこの表を取り込みます。
見てみましょう。
spline.c
/*
* spline.c:
* Spline interporator
* Coded by fanfan (S.Makino)
* Copyright (c) 2021 Agri-hitech LLC.
*/
#define REAL float
//#define REAL double
typedef struct {
REAL in;
REAL out;
REAL d;
} TBL;
#include "data.c"
/***********************************************************
spline.c -- スプライン補間
***********************************************************/
#define Nd sizeof(data)/sizeof(TBL)
float spline(float t, const TBL * T, int size)
{
int i, j, k;
REAL d, h;
i = 0; j = size - 1;
while (i < j) {
k = (i + j) / 2;
if (T[k].in < t) i = k + 1; else j = k;
}
if (i > 0) i--;
h = T[i + 1].in - T[i].in;
d = t - T[i].in;
return (((T[i + 1].d - T[i].d) * d / h + T[i].d * 3) * d
+ ((T[i + 1].out - T[i].out) / h
- (T[i].d * 2 + T[i + 1].d) * h)) * d + T[i].out;
}
#define DEBUG_MAIN
#ifdef DEBUG_MAIN
#include <stdio.h>
int main(){
int i;
REAL a;
for (i = 0; i <= 90; i++) {
a = (float) i / 10.0;
printf("%.16g,%.16g\n",a, spline(a, data, Nd));
}
}
#endif
#define REAL float と宣言してますが、RAM/ROMに余裕があるなら、#define REAL doubleに変えてもいいです。精度が上がりますが、普通のセンサーシステムでは不要な範囲です。
実行してみましょう。
gcc spline.c -o spline
./spline
結果をexcl等でグラフ化します。
まとめ
- 目標は達成されました
- 構造体の使い方の例になっていると思います
- patchをくださった方ありがとうございます。