Simulink ルックアップテーブルのCコード実装

1次元ルックアップテーブル

Simulink で生成される1次元ルックアップテーブルのコードは、以下のような構造体と関数で構成される。

typedef struct {
    double y_values[10]; // Y軸データ(出力値)
    double x_values[10]; // X軸データ(入力ブレークポイント)
} LookupTable1D_Params;

double interpolate_1d(double input, const double* x_bp, const double* y_table,
                      unsigned int max_idx) {
    double fraction;
    unsigned int left_idx;

    if (input <= x_bp[0]) {
        left_idx = 0;
        fraction = 0.0;
    } else if (input < x_bp[max_idx]) {
        unsigned int mid, right = max_idx, left = 0;
        mid = max_idx >> 1;
        while (right - left > 1) {
            if (input < x_bp[mid])
                right = mid;
            else
                left = mid;
            mid = (left + right) >> 1;
        }
        left_idx = left;
        fraction = (input - x_bp[left_idx]) / (x_bp[left_idx + 1] - x_bp[left_idx]);
    } else {
        left_idx = max_idx - 1;
        fraction = 1.0;
    }

    double y_low = y_table[left_idx];
    return y_low + (y_table[left_idx + 1] - y_low) * fraction;
}

この関数は、入力が範囲外の場合も補間を継続せず、端点での線形外挿(または固定値)ではなく、端点間の補間比率を0または1に固定することで、安全な挙動を保証する。

使用例:

#include <stdio.h>

const LookupTable1D_Params params = {
    .y_values = {1.2779, 2.40725, 2.65979, 2.23011, 1.58076,
                 1.20665, 0.98394, 0.91311, 1.01913, 1.14780},
    .x_values = {0.0, 0.33, 0.66, 0.99, 1.32, 1.65, 1.98, 2.31, 2.64, 2.97}
};

int main() {
    double out = interpolate_1d(0.25, params.x_values, params.y_values,
                                sizeof(params.x_values)/sizeof(params.x_values[0]) - 1);
    printf("出力: %f\n", out);
    return 0;
}

整数型対応版のカスタマイズ例

組込みシステム向けに、uint16_t を使用したバージョンも実装可能である。

#include <stdint.h>

uint16_t height_current_map[5][2] = {
    {100, 500},
    {500, 800},
    {2000, 1100},
    {3000, 1300},
    {4000, 1400}
};

uint16_t lookup_uint16(uint16_t in, uint16_t (*tbl)[2], unsigned int max_idx) {
    double frac;
    unsigned int left;

    if (in <= tbl[0][0]) {
        left = 0;
        frac = 0.0;
    } else if (in < tbl[max_idx][0]) {
        unsigned int mid, right = max_idx, l = 0;
        mid = max_idx >> 1;
        while (right - l > 1) {
            if (in < tbl[mid][0])
                right = mid;
            else
                l = mid;
            mid = (l + right) >> 1;
        }
        left = l;
        frac = (double)(in - tbl[left][0]) / (double)(tbl[left + 1][0] - tbl[left][0]);
    } else {
        left = max_idx - 1;
        frac = 1.0;
    }

    double y0 = tbl[left][1];
    return (uint16_t)(y0 + (tbl[left + 1][1] - y0) * frac);
}

// 使用例(省略)

2次元ルックアップテーブル

2次元テーブルでは、X軸とY軸のブレークポイントと、Z軸(出力)の格子データを持つ。メモリレイアウトは行優先(row-major)で、X軸が高速変化方向(stride = Xサイズ)となる。

typedef struct {
    double z_data[9];      // Z値(出力テーブル)
    double x_bp[3];        // X軸ブレークポイント
    double y_bp[3];        // Y軸ブレークポイント
    unsigned int max_idx[2]; // 各軸の最大インデックス(size - 1)
} LookupTable2D_Params;

double interpolate_2d(double x_in, double y_in,
                      const double* x_bp, const double* y_bp,
                      const double* z_table,
                      const unsigned int max_idx[2],
                      unsigned int x_stride) {
    double fx, fy;
    unsigned int ix, iy;

    // X方向の補間係数とインデックス
    if (x_in <= x_bp[0]) {
        ix = 0; fx = 0.0;
    } else if (x_in < x_bp[max_idx[0]]) {
        unsigned int l = 0, r = max_idx[0], m = max_idx[0] >> 1;
        while (r - l > 1) {
            if (x_in < x_bp[m]) r = m; else l = m;
            m = (l + r) >> 1;
        }
        ix = l;
        fx = (x_in - x_bp[ix]) / (x_bp[ix + 1] - x_bp[ix]);
    } else {
        ix = max_idx[0] - 1; fx = 1.0;
    }

    // Y方向の補間係数とインデックス
    if (y_in <= y_bp[0]) {
        iy = 0; fy = 0.0;
    } else if (y_in < y_bp[max_idx[1]]) {
        unsigned int l = 0, r = max_idx[1], m = max_idx[1] >> 1;
        while (r - l > 1) {
            if (y_in < y_bp[m]) r = m; else l = m;
            m = (l + r) >> 1;
        }
        iy = l;
        fy = (y_in - y_bp[iy]) / (y_bp[iy + 1] - y_bp[iy]);
    } else {
        iy = max_idx[1] - 1; fy = 1.0;
    }

    // 4点によるバイリニア補間
    size_t idx00 = iy * x_stride + ix;
    double z00 = z_table[idx00];
    double z10 = z_table[idx00 + 1];
    double z01 = z_table[idx00 + x_stride];
    double z11 = z_table[idx00 + x_stride + 1];

    double z_x0 = z00 + (z10 - z00) * fx;
    double z_x1 = z01 + (z11 - z01) * fx;

    return z_x0 + (z_x1 - z_x0) * fy;
}

使用例:

const LookupTable2D_Params tbl2d = {
    .z_data = {4.0, 16.0, 10.0, 5.0, 19.0, 18.0, 6.0, 20.0, 23.0},
    .x_bp = {1.0, 2.0, 3.0},
    .y_bp = {4.0, 5.0, 6.0},
    .max_idx = {2, 2}
};

int main() {
    double result = interpolate_2d(1.5, 5.1,
                                   tbl2d.x_bp, tbl2d.y_bp, tbl2d.z_data,
                                   tbl2d.max_idx, 3);
    printf("2D出力: %f\n", result);
    return 0;
}

非正方グリッドにも対応可能であり、x_stride をX軸の要素数に設定することで、任意の矩形テーブルに対応できる。たとえばX軸が4点、Y軸が3点の場合、max_idx = {3, 2}x_stride = 4 とする。

タグ: Simulink C Embedded Systems Lookup Table Interpolation

8月8日 16:01 投稿