こんにちは。
仕事の自動化にやりがいと達成感を感じるガッくんです。

この記事の目次

背景・目的

高速フーリエ変換（ FFT ）をやっていて精度がイマイチだったので、色々考えているときに、波形を近似してデータ数を増やすと精度が上がったりしないかなと何の根拠もなしに思いました。

FFT の精度が上がるかどうかはとりあえず置いといて、補間方法をまとめることにしました。

動作環境

・Windows 7
・winpython 64bit 3.4.4

プログラム

ソースコード


###############################################################################
# 曲線を近似するプログラム
###############################################################################
#==============================================================================
# ライブラリインポート
#==============================================================================
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
 
#==============================================================================
# 1 次スプライン補間関数
#==============================================================================
def spline_primary(x, y, point):
    f = interpolate.interp1d(x, y)                         # 補間
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    Y = f(X)                                               # 準備した x に対応する y の計算
    return X, Y                                            # 戻り値
 
#==============================================================================
# 2 次スプライン補間関数
#==============================================================================
def spline_secondary(x, y, point):
    f = interpolate.interp1d(x, y, kind="quadratic")       # 補間
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    Y = f(X)                                               # 準備した x に対応する y の計算
    return X, Y                                            # 戻り値
 
#==============================================================================
# 3 次スプライン補間関数
#==============================================================================
def spline_3rd_order(x, y, point):
    f = interpolate.interp1d(x, y, kind="cubic")           # 補間
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    Y = f(X)                                               # 準備した x に対応する y の計算
    return X, Y                                            # 戻り値
 
#==============================================================================
# 秋間補間関数
#==============================================================================
def spline_akima(x, y, point):
    f = interpolate.Akima1DInterpolator(x, y)              # 補間
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    Y = f(X)                                               # 準備した x に対応する y の計算
    return X, Y                                            # 戻り値
 
#==============================================================================
# グラフ描画関数
#==============================================================================
def graph(x, y, a, b, interpolation_method, xmin, xmax, ymin, ymax, x_scale_interval, y_scale_interval):
    # 元々の値の出力
    plt.plot(x, y, 'ro', label="data point")
   
    # 補間した結果を出力
    plt.plot(a, b, label=interpolation_method)
   
    # グラフのタイトル・目盛設定
    plt.title("spline")
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
   
    # 凡例の位置設定
    plt.legend(loc='lower right')
   
    # 格子の設定
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
    plt.xticks(list(filter(lambda x: x%x_scale_interval==0, np.arange(xmin, xmax))))
    plt.yticks(list(filter(lambda x: x%y_scale_interval==0, np.arange(ymin, ymax))))
   
    # グラフ出力
    file_name = interpolation_method + '.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
 
#==============================================================================
# テスト
#==============================================================================
# 値の準備
x = [-9, -5, -2,  1,  3,  7, 12, 16]
y = [-1, -4,  2, -2, -4,  0,  4, 10]
 
# 補間
interpolation_method1 = 'spline 1st order'
a1, b1 = spline_primary(x, y, 100)   # 1 次スプライン補間
 
interpolation_method2 = 'spline 2nd order'
a2, b2 = spline_secondary(x, y, 100) # 2 次スプライン補間
 
interpolation_method3 = 'spline 3rd order'
a3, b3 = spline_3rd_order(x, y, 100) # 3 次スプライン補間
 
interpolation_methoda = 'akima'
aa, ba = spline_akima(x, y, 100)     # 秋間スプライン補間
 
# グラフ出力設定
xmin = min(x) - 1
xmax = max(x) + 1
ymin = min(y) - 5
ymax = max(y) + 1
x_scale_interval = 2 # x 方向目盛間隔
y_scale_interval = 2 # y 方向目盛間隔
 
# グラフ出力（ 1 ケースずつ）
graph(x, y, a1, b1, interpolation_method1, xmin, xmax, ymin, ymax, x_scale_interval, y_scale_interval)
graph(x, y, a2, b2, interpolation_method2, xmin, xmax, ymin, ymax, x_scale_interval, y_scale_interval)
graph(x, y, a3, b3, interpolation_method3, xmin, xmax, ymin, ymax, x_scale_interval, y_scale_interval)
graph(x, y, aa, ba, interpolation_methoda, xmin, xmax, ymin, ymax, x_scale_interval, y_scale_interval)
 
# グラフ出力（ 2 次、 3 次、秋間の重ね合わせ） ========================================
# 元々の値の出力
plt.plot(x, y, 'ro', label="data point")
 
# 補間した結果を出力
plt.plot(a2, b2, label=interpolation_method2)
plt.plot(a3, b3, label=interpolation_method3)
plt.plot(aa, ba, label=interpolation_methoda)
 
# グラフのタイトル・目盛設定
plt.title("spline")
plt.xlim([xmin, xmax])
plt.ylim([ymin, ymax])
 
# 凡例の位置設定
plt.legend(loc='lower right')
 
# 格子の設定
plt.grid(which='major', color='black', linestyle='-')
plt.grid(which='minor', color='black', linestyle='-')
plt.xticks(list(filter(lambda x: x%x_scale_interval==0, np.arange(xmin, xmax))))
plt.yticks(list(filter(lambda x: x%y_scale_interval==0, np.arange(ymin, ymax))))
 
# グラフ出力
file_name = 'comparison interpolation method.jpg'
plt.savefig(file_name)
 
# 表示
plt.show()

結果

高精度な補間を探していて、波形の様なデータにはスプライン補間が良さそうだったので、 1 次〜 3 次スプライン補間と秋間補間を試しました。

1 次スプラインは直線でデータを結んでいるだけなので、使うことはないと思うので放置します。

2 次スプラインと 3 次スプラインでは 3 次スプラインの方がほんの少し滑らかになっています。
秋間補間は 2 次スプラインと 3 次スプラインよりもさらに滑らかになっています。

基本的にスプライン補間はデータ全点を通るのですが、無理矢理全点を通る曲線にするため、振動しやすいみたいです。

秋間補間は振動が極端に小さいのが特徴ですね。