cvl-robot's diary

研究ノート メモメモ https://github.com/dotchang/

トラ技2018年10月号で紹介されている無限スプライン関数の計算アルゴリズムを試してみる

データストリームに対してリアルタイムで逐次に3次スプライン補間を計算できるアルゴリズムだそうです。

トランジスタ技術 2018年 10 月号

トランジスタ技術 2018年 10 月号

SSDAC特設サイト

#include <iostream>
#include <math.h>
#include <vector>
#include <iterator>
#include <algorithm>

#include <omp.h>

double calc_coefficient_b(double *d, int n)
{
  static const double alpha = -2. + sqrt(3.);
  static const double a = -3. * (sqrt(3.) - 1.);
  static const double b = 3. * (2. * sqrt(3.) - 3.);
  double e(0.), f(0.);
  for(int i=n; i>0; i--){
    f = d[n+i] + alpha * f;
    e = d[n-i] + alpha * e;
  }
  return a * d[n] + b * (f + e);
}

int main() {
  static const int n = 9;
  static const int multiples = 64;
  std::vector<double> output;
  output.clear();

  // input-data preparation
  double samples[] = {0., 2., 3., 5., 0., -1., 2., -3., -4., -3., 0., 3., 3., 6., 0., 9., 1.};
  std::vector<double> data(samples, std::end(samples));
  int n_samples = data.size();
  data.insert(data.begin(), n, 0.);
  data.insert(data.end(), n, 0.);

  // calc-preparation
  std::vector<double> x3, x2, x1;
  x3.clear();
  x2.clear();
  x1.clear();
  for(int i = 0; i < multiples; i++){
    double x = static_cast<double>(i) / static_cast<double>(multiples);
    x1.push_back(x);
    x2.push_back(pow(x, 2));
    x3.push_back(pow(x, 3));
  }
  std::vector<double> interpolated;
  interpolated.resize(multiples);
  
  int start_idx = 0;
  double bj, bj_1, dj, dj_1;
  bj = calc_coefficient_b(&data[start_idx], n);
  dj = data[start_idx];
  for(int j = start_idx; j < n_samples; j++){
    // calc coefficients
    bj_1 = calc_coefficient_b(&data[j + 1], n);
    dj_1 = data[j + n + 1];
    double aj = (bj_1 - bj) / 3.;
    double cj = dj_1 - dj - 2. * bj / 3. - bj_1 / 3.;

    // calc oversampling
#pragma omp parallel for
    for(int i = 0; i < multiples; i++){
      interpolated[i] = aj * x3[i] + bj * x2[i] + cj * x1[i] + dj;
    }
    std::copy(interpolated.begin(), interpolated.end(), std::back_inserter(output));
    
    bj = bj_1;
    dj = dj_1;
  }

  // outputs
  std::cout << "samples:" << std::endl;
  for(int i = start_idx; i < n_samples; i++){
    std::cout << samples[i] << std::endl;
    for(int j = 0; j < multiples - 1; j++){
      std::cout << 0. << std::endl;
    }
  }
  std::cout << "output:" << std::endl;
  for(int i=0; i<output.size(); i++){
    std::cout << output[i] << std::endl;
  }

  return 0;
}

適当な入力データに対する計算結果はこんな感じ。

b:
2.1531

  • 2.03546

2.98887

  • 6.91988

3.69088
4.1564

  • 8.31614

5.10831

  • 0.117074

1.35999
0.677115

  • 4.06843

6.59671

  • 13.3183

19.6765

  • 20.3877

10.8744

  • 2.10998

f:id:cvl-robot:20180910174733p:plain

きれいにスプライン補間できていそうですね。
あまりにもアルゴリズムが簡単すぎて、これが本当に新しくて凄いものなのかピンと来ていません。

[1] https://www.tezukuri-amp.org/bunkakai/dac/cgi-bin/img-box/img20170804205124.pdf

GitHub - kritzikratzi/ofxAvCodec: openFrameworks wrapper for FFmpeg/libavcodec
GitHub - roymacdonald/ofxSoundObjects: openFrameworks addon implementation for the ofSoundObject implementaion developed at the YCAM '13 ofDevCon, yet never released