解析エンジニアの自動化 blog

コツコツと自動化した方法を残す blog

Python で曲面近似(サーフェスフィッティング)する



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


この記事の目次



背景・目的


曲線近似はわりと簡単に出来ました。

もちろん、簡単な近似プログラムなので出来たというのは大げさかもしれませんが…

曲線近似をやったので曲面近似(サーフェスフィッティング)もやってみたいと思います。

リンクの記事の曲線近似と計算手法自体は同じです。
Python で 放射基底関数 (RBF) 補間 - 解析エンジニアの自動化 blog

テスト用のある程度複雑な曲面を作り、 その曲面上の点をランダムに抽出します。
今回は誤差は与えません。

そして、 ランダムに抽出した曲面上の点に対して近似曲面を作ります。



動作環境


Windows 10
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 曲面近似プログラム
###############################################################################
#==============================================================================
# ライブラリインポート
#==============================================================================
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
#==============================================================================
# 曲面作成
#==============================================================================
x, y, z = 10 * np.random.random((3, 10)) # 乱数で適当に x, y, z を作成する
xi, yi = np.linspace(x.min(), x.max(), 100), np.linspace(y.min(), y.max(), 100) # x, y の値の最小・最大から等間隔の配列を作成する
xi, yi = np.meshgrid(xi, yi) # xi, yi から mesh を作成する
 
# カーネル種類
#   線形 RBF => linear
#   ガウシアン RBF => gaussian
#   多重二乗 RBF => multiquadric
#   逆二乗 RBF => inverse
#   多重調和スプライン RBF => cubic(3次), quintic(5次), thin_plate(薄板スプライン)
rbf = interpolate.Rbf(x, y, z, function='gaussian') # x, y, z の値で RBF 補間をして曲面を作成する
 
zi = rbf(xi, yi) # x, y, z の曲面における xi, yi の位置の zi を計算する
 
#==============================================================================
# 3 次元グラフ
#==============================================================================
fig = plt.figure(figsize=(18, 7), dpi=200) # 画像を作成する
el = 55                                    # 視点高さを設定する
 
#==============================================================================
ax = fig.add_subplot(231, projection='3d') # 2 × 3 の 1 枚目に描画する
 
# 曲面のプロット
ax.plot_surface(xi, yi, zi)                      # サーフェスの描画
ax.view_init(elev=el, azim=10)                   # ビューの設定
ax.set_title('elev = ' + str(el) + ', deg = 10') # タイトルの設定
ax.set_xlabel('xi')                              # 軸ラベルの設定
ax.set_ylabel('yi')                              # 軸ラベルの設定
ax.set_zlabel('zi')                              # 軸ラベルの設定
 
#==============================================================================
ax = fig.add_subplot(232, projection='3d') # 2 × 3 の 2 枚目に描画する
 
# 曲面のプロット
ax.plot_surface(xi, yi, zi)                      # サーフェスの描画
ax.view_init(elev=el, azim=45)                   # ビューの設定
ax.set_title('elev = ' + str(el) + ', deg = 45') # タイトルの設定
ax.set_xlabel('xi')                              # 軸ラベルの設定
ax.set_ylabel('yi')                              # 軸ラベルの設定
ax.set_zlabel('zi')                              # 軸ラベルの設定
 
#==============================================================================
ax = fig.add_subplot(233, projection='3d') # 2 × 3 の 3 枚目に描画する
 
# 曲面のプロット
ax.plot_surface(xi, yi, zi)                      # サーフェスの描画
ax.view_init(elev=el, azim=80)                   # ビューの設定
ax.set_title('elev = ' + str(el) + ', deg = 80') # タイトルの設定
ax.set_xlabel('xi')                              # 軸ラベルの設定
ax.set_ylabel('yi')                              # 軸ラベルの設定
ax.set_zlabel('zi')                              # 軸ラベルの設定
 
#==============================================================================
ax = fig.add_subplot(234, projection='3d') # 2 × 3 の 4 枚目に描画する
 
# 曲面のプロット
#   rstride と cstride はステップサイズ
#   cmap は彩色
#   linewidth は曲面のメッシュの線の太さ
ax.plot_wireframe(xi, yi, zi, rstride=1, cstride=1, cmap='hsv', linewidth=0.2) # ワイヤーフレームの描画
ax.view_init(elev=el, azim=10) # ビューの設定
ax.set_xlabel('xi')            # 軸ラベルの設定
ax.set_ylabel('yi')            # 軸ラベルの設定
ax.set_zlabel('zi')            # 軸ラベルの設定
 
#==============================================================================
ax = fig.add_subplot(235, projection='3d') # 2 × 3 の 5 枚目に描画する
 
# 曲面のプロット
#   rstride と cstride はステップサイズ
#   cmap は彩色
#   linewidth は曲面のメッシュの線の太さ
ax.plot_wireframe(xi, yi, zi, rstride=1, cstride=1, cmap='hsv', linewidth=0.2) # ワイヤーフレームの描画
ax.view_init(elev=el, azim=45) # ビューの設定
ax.set_xlabel('xi')            # 軸ラベルの設定
ax.set_ylabel('yi')            # 軸ラベルの設定
ax.set_zlabel('zi')            # 軸ラベルの設定
 
#==============================================================================
ax = fig.add_subplot(236, projection='3d') # 2 × 3 の 6 枚目に描画する
 
# 曲面のプロット
#   rstride と cstride はステップサイズ
#   cmap は彩色
#   linewidth は曲面のメッシュの線の太さ
ax.plot_wireframe(xi, yi, zi, rstride=1, cstride=1, cmap='hsv', linewidth=0.2) # ワイヤーフレームの描画
ax.view_init(elev=el, azim=80) # ビューの設定
ax.set_xlabel('xi')            # 軸ラベルの設定
ax.set_ylabel('yi')            # 軸ラベルの設定
ax.set_zlabel('zi')            # 軸ラベルの設定
 
#==============================================================================
# グラフ出力
file_name = 'Various 3D Images.jpg' # グラフ名設定
plt.savefig(file_name)              # グラフ出力
 
plt.show() # 描画
 
#==============================================================================
# 2 次元グラフ
#==============================================================================
fig = plt.figure(figsize=(18, 7), dpi=200) # 画像を作成する
 
# 曲面のプロット ==================================================================
ax = fig.add_subplot(121, projection='3d') # 1 × 2 の 1 枚目に描画する
ax.plot_surface(xi, yi, zi)                # サーフェスの描画
el = 100                                   # 視点高さを設定する
ax.view_init(elev=el, azim=90)             # ビューの設定
ax.set_title('elev = 100, deg = 90')       # タイトルの設定
ax.set_xlabel('xi')                        # 軸ラベルの設定
ax.set_ylabel('yi')                        # 軸ラベルの設定
ax.set_zlabel('zi')                        # 軸ラベルの設定
 
# コンター =======================================================================
ax = fig.add_subplot(122) # 1 × 2 の 2 枚目に描画する
contour = ax.contourf(xi, yi, zi)
fig.colorbar(contour)
ax.set_xlim([x.max(), x.min()])
ax.set_ylim([y.max(), y.min()])
ax.set_title('contour')       # タイトルの設定
ax.set_xlabel('xi')            # 軸ラベルの設定
ax.set_ylabel('yi')            # 軸ラベルの設定
 
# グラフ出力
file_name = '3D and Contour.jpg' # グラフ名設定
plt.savefig(file_name)           # グラフ出力
 
plt.show()



結果

このプログラムから出力される近似曲線の 3 次元グラフを図1 〜 図2 に示します。

図1 グラフ


図2 コンター




コメント

こんなに綺麗な画像が簡単に作れるとは思いませんでした。

曲線近似は使い所がマニアックなものしか思いつかないので、あまり今後はブログに出てこないかもしれません。



以上

Python で波形データを微分する



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


この記事の目次



背景・目的


積分をやったんで、微分も記事にまとめます。

2波形について試してみました。



動作環境


Windows 10
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 離散データの微分を計算するプログラム
###############################################################################
#==============================================================================
# ライブラリインポート
#==============================================================================
import numpy as np
import matplotlib.pyplot as plt
 
#==============================================================================
# テスト用離散データの作成
#==============================================================================
flag = 1
 
if flag == 0:
    x = np.linspace(0, 2.2, num=2**5+1)
    y = 2.1*x**4 - 3.3*x**3 - 2*x + 1
elif flag == 1:
    x = np.arange(-7.0, 8.0, 1)
    y = x**2
 
#==============================================================================
# 微分計算
#==============================================================================
dy = np.gradient(y)
 
#==============================================================================
# グラフ
#==============================================================================
# 微分対象波形
plt.plot(x, y, 'r-', label='wave')
plt.plot(x, y, 'r.')
 
# 微分波形
plt.plot(x, dy, 'b-', label='differential')
plt.plot(x, dy, 'b.')
 
# グラフのタイトル・目盛設定
plt.title("Differential")
 
# 凡例の位置設定
plt.legend(loc='upper left')
 
# 罫線
plt.grid(which='major', color='black', linestyle='-')
plt.grid(which='minor', color='black', linestyle='-')
 
# 軸ラベル
plt.xlabel('x')
plt.ylabel('y')
 
# y 軸目盛
if flag == 1:
    plt.ylim([-20, 70])
 
file_name = 'differential.jpg'
plt.savefig(file_name)
plt.show()



結果

ソースコードの『テスト用離散データの作成』の flag == 0 の時の微分波形を図1 、 flag == 1 の時の微分波形を図2 に示します。

図1 flag == 0


図2 flag == 1



コメント

図1 では微分が正しく行われているか分からなかったので、図2 から微分が正しく行われていることを確認しました。

これで物理量の変換も出来るようになりました。



以上

Python で波形データを積分する



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


この記事の目次



背景・目的


今まで波形の変換、近似、補間を記事にまとめてきました。

今度は変換とはまた違う " 積分 " を行なってみます。

積分すると物理量が変わり、波形の意味が変わります。

1つの波形から色々な情報を読み取るためにも積分についてまとめます。



動作環境


Windows 10
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 離散データの積分値の累積を計算するプログラム
###############################################################################
#==============================================================================
# ライブラリインポート
#==============================================================================
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
 
#==============================================================================
# テスト用離散データの作成
#==============================================================================
x = np.linspace(0, 2.2, num=2**5+1)
y = 2.1*x**4 - 3.3*x**3 - 2*x + 1
 
#==============================================================================
# 積分計算
#==============================================================================
integral_y_cumtrapz = integrate.cumtrapz(y, x, initial=0)
 
#==============================================================================
# グラフ
#==============================================================================
# 積分対象波形
plt.plot(x, y, 'r-', label='wave')
plt.plot(x, y, 'r.')
 
# 積分波形
plt.plot(x, integral_y_cumtrapz, 'b-', label='integral')
plt.plot(x, integral_y_cumtrapz, 'b.')
 
# グラフのタイトル・目盛設定
plt.title("Integral")
 
# 凡例の位置設定
plt.legend(loc='upper left')
 
# 罫線
plt.grid(which='major', color='black', linestyle='-')
plt.grid(which='minor', color='black', linestyle='-')
 
# 軸ラベル
plt.xlabel('x')
plt.ylabel('y')
 
file_name = 'integral.jpg'
plt.savefig(file_name)
plt.show()



結果

積分した結果を図1 に示します。

図1 積分


赤い波形を積分すると青い波形になりました。



コメント

波形データって奥が深いですね。

微分も後でやってみます。



以上

Python でチェビシェフ(Chebyshev)基底を使った高次の多項式近似



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


この記事の目次



背景・目的


前回の記事で polyfit で多項式近似をやりました。

Python で多項式近似と決定係数 - 解析エンジニアの自動化 blog

近似曲線を作る際になかなか便利なものができたと思います。
しかし、高次の多項式近似で Warning が出ていました。。

調べていくとチェビシェフ(Chebyshev)基底を使った多項式近似という方法が高次の曲線でも安定的に近似曲線が得られるらしいことが分かりました。

テスト用にある程度複雑な波形を作り、 誤差を与えます。

そして、誤差を含む波形上の点に対して近似曲線を作ります。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# チェビシェフ基底多項式近似
###############################################################################
#==============================================================================
# ファイルインポート
#==============================================================================
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Chebyshev as c
 
#==============================================================================
# グラフ描画関数
#==============================================================================
def wave_graph(x, y, x1, y1, xmin, xmax, ymin, ymax):
    plt.figure(figsize=(8, 6))
   
    # 元々の値の出力
    plt.plot(x, y, 'r.', label="original data point")
    plt.plot(x, y, 'r')
   
    # 誤差を含んだ値の出力
    plt.plot(x1, y1, 'b.', label="data point for test")
    plt.plot(x1, y1, 'b')
   
    # グラフのタイトル・目盛設定
    plt.title("wave for test")
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
   
    # 凡例の位置設定
    plt.legend(loc='upper right')
   
    # 格子の設定
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
   
    # グラフ出力
    file_name = 'wave for test.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
   
    return file_name
 
def err_graph(x, y, xmin, xmax, ymin, ymax):
    plt.figure(figsize=(8, 6))
   
    # 各 x 値における誤差の出力
    plt.bar(x, y, linewidth = 0, width = 0.0007, color="#1E7F00")
    plt.plot(x, y, 'r.', label="error of data point")
    plt.plot(x, y, 'r')
   
    # グラフのタイトル・目盛設定
    plt.title("error of wave for test")
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
   
    # 凡例の位置設定
    plt.legend(loc='upper right')
   
    # 誤差の表示
    max_err = max(y)
    min_err = min(y)
    plt.text(xmin+np.abs((xmax-xmin)*0.01), ymin+np.abs((ymax-ymin)*0.07), 'upper error = ' + str(max_err))
    plt.text(xmin+np.abs((xmax-xmin)*0.01), ymin+np.abs((ymax-ymin)*0.01), 'lower error = ' + str(min_err))
   
    # 格子の設定
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
   
    # グラフ出力
    file_name = 'error of wave for test.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
   
    return file_name
 
def graph(x, y, a, b, interpolation_method):
    plt.figure(figsize=(8, 6))
   
    # 元々の値の出力
    plt.plot(x, y, 'ro', label="data point")
   
    # 補間した結果を出力
    plt.plot(a, b, 'k', label=interpolation_method)
   
    # グラフのタイトル・目盛設定
    plt.title(interpolation_method)
    plt.grid(True)
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # 凡例の位置設定
    plt.legend(loc='lower right')
   
    # グラフ出力
    file_name = interpolation_method + '.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
 
#==============================================================================
# チェビシェフ基底を使った多項式近似関数
#==============================================================================
def chebyshev(x, y, stop_n, point):
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    max_corrcoef = 0                                       # 最終的に採用した決定係数記憶用変数
    orders = []                                            # 次数の履歴記憶用変数
    corrcoefs = []                                         # 決定係数の履歴記憶用変数
    eps = 0.00001                                          # 打切り誤差
    for n in range(1, stop_n):                             # 1 ~ stop_n まで繰り返す
        f = c.fit(x, y, n)                                 # 多項式近似
        Y = f(X)                                           # X に対応する y の値を計算する
        corrcoef = np.corrcoef(y, f(x))[0, 1]**2           # 決定係数を計算する
        orders.append(n)                                   # 次数の履歴記憶
        corrcoefs.append(corrcoef)                         # 決定係数の履歴記憶
        # 現状最大の決定係数 max_corrcoef とさっき計算した corrcoef の差の絶対値が打切り誤差より小さい場合
        if abs(max_corrcoef - corrcoef) < eps:
            break                                          # 多項式計算のループ終了
        # 現状最大の決定係数 max_corrcoef よりさっき計算した corrcoef が大きい場合
        if max_corrcoef < corrcoef:
            max_corrcoef = corrcoef                        # max_corrcoef の更新
            fix_Y = Y                                      # 最大決定係数の時の Y を記憶する
            fix_n = n                                      # 最大決定係数を記憶する
    plt.figure(figsize=(8, 6))                             # グラフの初期化
    plt.plot(orders, corrcoefs, 'r.', label="coefficient of determination") # 決定係数の値の出力
    plt.plot(orders, corrcoefs, 'r')                       # 値のドットを線で結ぶ
    plt.legend(loc='lower right')                          # 凡例の位置設定
    plt.grid(True)                                         # 罫線表示
    plt.xlabel('order')                                    # 横軸タイトル
    plt.ylabel('coefficient of determination')             # 縦軸タイトル
    file_name = 'Determination coefficient transition to ' + str(n) + '.jpg' # グラフ名設定
    plt.savefig(file_name)                                 # グラフ出力
    return X, fix_Y, fix_n
 
#==============================================================================
# 決定係数が最大を示す n 次式の最小二乗法近似関数
#==============================================================================
def n_th_order_least_squares_approximation_by_max_corrcoef(x, y, stop_n, point):
   
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    max_corrcoef = 0                                       # 最終的に採用した決定係数記憶用変数
    orders = []                                            # 次数の履歴記憶用変数
    corrcoefs = []                                         # 決定係数の履歴記憶用変数
    for n in range(1, stop_n):                             # 1 ~ stop_n まで繰り返す
        f = np.poly1d(np.polyfit(x, y, n))                 # 多項式近似
        Y = f(X)                                           # X に対応する y の値を計算する
        corrcoef = np.corrcoef(y, f(x))[0, 1]**2           # 決定係数を計算する
        orders.append(n)                                   # 次数の履歴記憶
        corrcoefs.append(corrcoef)                         # 決定係数の履歴記憶
        # 現状最大の決定係数 max_corrcoef よりさっき計算した corrcoef が大きい場合
        if max_corrcoef < corrcoef:
            max_corrcoef = corrcoef                        # max_corrcoef の更新
            fix_Y = Y                                      # 最大決定係数の時の Y を記憶する
            fix_n = n                                      # 最大決定係数を記憶する
    plt.figure(figsize=(8, 6))                             # グラフの初期化
    plt.plot(orders, corrcoefs, 'r.', label="coefficient of determination") # 決定係数の値の出力
    plt.plot(orders, corrcoefs, 'r')                       # 値のドットを線で結ぶ
    plt.legend(loc='lower right')                          # 凡例の位置設定
    plt.grid(True)                                         # 罫線表示
    plt.xlabel('order')                                    # 横軸タイトル
    plt.ylabel('coefficient of determination')             # 縦軸タイトル
    file_name = 'Determination coefficient transition to ' + str(fix_n) + '.jpg' # グラフ名設定
    plt.savefig(file_name)                                 # グラフ出力
    return X, fix_Y, fix_n                                 # 戻り値
 
#==============================================================================
# 決定係数が上昇しなくなった n 次式の最小二乗法近似関数
#==============================================================================
def n_th_order_least_squares_approximation_by_eps(x, y, stop_n, point):
   
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    max_corrcoef = 0                                       # 最終的に採用した決定係数記憶用変数
    orders = []                                            # 次数の履歴記憶用変数
    corrcoefs = []                                         # 決定係数の履歴記憶用変数
    eps = 0.00001                                          # 打切り誤差
    for n in range(1, stop_n):                             # 1 ~ stop_n まで繰り返す
        f = np.poly1d(np.polyfit(x, y, n))                 # 多項式近似
        Y = f(X)                                           # X に対応する y の値を計算する
        corrcoef = np.corrcoef(y, f(x))[0, 1]**2           # 決定係数を計算する
        orders.append(n)                                   # 次数の履歴記憶
        corrcoefs.append(corrcoef)                         # 決定係数の履歴記憶
        # 現状最大の決定係数 max_corrcoef とさっき計算した corrcoef の差の絶対値が打切り誤差より小さい場合
        if abs(max_corrcoef - corrcoef) < eps:
            break                                          # 多項式計算のループ終了
        # 現状最大の決定係数 max_corrcoef よりさっき計算した corrcoef が大きい場合
        if max_corrcoef < corrcoef:
            max_corrcoef = corrcoef                        # max_corrcoef の更新
            fix_Y = Y                                      # 最大決定係数の時の Y を記憶する
            fix_n = n                                      # 最大決定係数を記憶する
    plt.figure(figsize=(8, 6))                             # グラフの初期化
    plt.plot(orders, corrcoefs, 'r.', label="coefficient of determination") # 決定係数の値の出力
    plt.plot(orders, corrcoefs, 'r')                       # 値のドットを線で結ぶ
    plt.legend(loc='lower right')                          # 凡例の位置設定
    plt.grid(True)                                         # 罫線表示
    plt.xlabel('order')                                    # 横軸タイトル
    plt.ylabel('coefficient of determination')             # 縦軸タイトル
    file_name = 'Determination coefficient transition to ' + str(fix_n) + '.jpg' # グラフ名設定
    plt.savefig(file_name)                                 # グラフ出力
    return X, fix_Y, fix_n                                 # 戻り値
 
#==============================================================================
# テスト用データ作成
#==============================================================================
# 値の準備 =====================================================================
N = 2**8                              # サンプリング数
dt = 0.0009                           # サンプリング周期
ampl = np.array([   1,   1,   2,  1]) # 振幅
freq = np.array([  10,  20,  30, 40]) # 周波数
peri_axis = np.arange(0, N*dt, dt)    # 周期軸
error = 0.11                          # 誤差係数
for i in range(len(freq)):
    if i == 0:
        origin_y = ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis)                       # 入力波形
        y = ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis) + error * np.random.randn(N) # 誤差を与えた入力波形
    else:
        origin_y += ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis)                       # 入力波形
        y += ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis) + error * np.random.randn(N) # 誤差を与えた入力波形
 
# 与えた誤差の計算 ==============================================================
err = (y) - (origin_y)
 
# グラフ出力設定 ================================================================
xmin = min(peri_axis) - (3*dt) # グラフの x 軸の最小値設定
xmax = max(peri_axis) + (3*dt) # グラフの x 軸の最大値設定
ymin = min(y) * 1.2            # グラフの y 軸の最小値設定
ymax = max(y) * 1.2            # グラフの y 軸の最大値設定
x_scale_interval = 1           # x 方向目盛間隔
y_scale_interval = 1           # y 方向目盛間隔
 
max_err = max(err)
min_err = min(err)
ymin_err = min_err * 1.5
ymax_err = max_err * 1.5
 
# グラフ出力 ====================================================================
wave_graph(peri_axis, origin_y, peri_axis, y, xmin, xmax, ymin, ymax)
err_graph(peri_axis, err, xmin, xmax, ymin_err, ymax_err)
 
# 決定係数が最大を示す次数による多項式補間 =========================================
X, Y, n = n_th_order_least_squares_approximation_by_max_corrcoef(peri_axis, y, 200, 1000)
graph(peri_axis, origin_y, X, Y, str(n) + 'th order by max coef')
 
# 決定係数が上昇しなくなった次数による多項式補間 ======================================
X, Y, n = n_th_order_least_squares_approximation_by_eps(peri_axis, y, 200, 1000)
graph(peri_axis, origin_y, X, Y, str(n) + 'th order by stop coef')
 
# チェビシェフ基底を使った多項式近似 =================================================
X, Y, n = chebyshev(peri_axis, y, 200, 1000)
graph(peri_axis, origin_y, X, Y, 'chebyshev by ' + str(n) + 'th order')



結果

polyfit による多項式近似曲線は前回の記事で作成済みなので、流用して、チェビシェフ基底を使った多項式近似曲線と比較してみました。

具体的には下記の①〜③の近似曲線を作成して傾向を確認しました。

①決定係数が最大を示した次数の近似曲線
②決定係数の上昇が無くなった次数の近似曲線
③ 定係数の上昇が無くなった次数のチェビシェフ近似曲線

図1 作成波形グラフ


図2 誤差グラフ


図3 ①の決定係数グラフ


図4 ①の近似曲線グラフ


図5 ②の決定係数グラフ


図6 ②の近似曲線グラフ


図7 ③の決定係数グラフ


図8 ③の近似曲線グラフ


①の近似曲線は 194 次という大きな次数の近似曲線です。
図3 から次数を上げるにつれて決定係数も上昇している傾向が分かります。
元々の波形が複雑なので決定係数が 1 付近に来る決定係数の傾向も乱れが生じています。

図4 は 0.24 秒あたりで波形が大きく乱れています。
決定係数は高いのですが、全く近似出来ていません。

図5 を見ると決定係数がほぼ一定になる領域に達してすぐに終了しています。

図6 は図4 程は乱れが無いですが、近似精度は良くありません。

やはり、決定係数の値が 1 に近ければ良いのではなく、決定係数が上昇傾向から一定に移り変わる時の次数による近似曲線のほうが良く曲線を表現出来そうです。

図7 を見ると図3 、図5 とは違い、 1 に近い一定となる領域での乱れがありません。

図8 の近似曲線はかなり精度よく曲線を表現出来ています。



コメント

チェビシェフ基底を使った多項式近似では Warning はスローされませんでした。

多項式系の近似曲線は精度さえ良ければ大幅に欠損したデータの補間なども使えそうです。



以上

Python で多項式近似と決定係数



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


この記事の目次



背景・目的


今まで近似や補間をブログに書いてきましたが、多項式近似を忘れていました。
1 番ポピュラーな近似方法を忘れていました。

多項式の近似は案外使い勝手の良いシロモノです。

テスト用にある程度複雑な波形を作り、 誤差を与えます。

そして、誤差を含む波形上の点に対して近似曲線を作ります。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# n 次多項式補間
###############################################################################
#==============================================================================
# ファイルインポート
#==============================================================================
import numpy as np
import matplotlib.pyplot as plt
 
#==============================================================================
# グラフ描画関数
#==============================================================================
def wave_graph(x, y, x1, y1, xmin, xmax, ymin, ymax):
    plt.figure(figsize=(8, 6))
   
    # 元々の値の出力
    plt.plot(x, y, 'r.', label="original data point")
    plt.plot(x, y, 'r')
   
    # 誤差を含んだ値の出力
    plt.plot(x1, y1, 'b.', label="data point for test")
    plt.plot(x1, y1, 'b')
   
    # グラフのタイトル・目盛設定
    plt.title("wave for test")
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
   
    # 凡例の位置設定
    plt.legend(loc='upper right')
   
    # 格子の設定
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
   
    # グラフ出力
    file_name = 'wave for test.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
   
    return file_name
 
def err_graph(x, y, xmin, xmax, ymin, ymax):
    plt.figure(figsize=(8, 6))
   
    # 各 x 値における誤差の出力
    plt.bar(x, y, linewidth = 0, width = 0.0007)
    plt.plot(x, y, 'r.', label="error of data point")
    plt.plot(x, y, 'r')
   
    # グラフのタイトル・目盛設定
    plt.title("error of wave for test")
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
   
    # 凡例の位置設定
    plt.legend(loc='upper right')
   
    # 誤差の表示
    max_err = max(y)
    min_err = min(y)
    plt.text(xmin+np.abs((xmax-xmin)*0.01), ymin+np.abs((ymax-ymin)*0.07), 'upper error = ' + str(max_err))
    plt.text(xmin+np.abs((xmax-xmin)*0.01), ymin+np.abs((ymax-ymin)*0.01), 'lower error = ' + str(min_err))
   
    # 格子の設定
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
   
    # グラフ出力
    file_name = 'error of wave for test.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
   
    return file_name
 
def graph(x, y, a, b, interpolation_method):
    plt.figure(figsize=(8, 6))
   
    # 元々の値の出力
    plt.plot(x, y, 'ro', label="data point")
   
    # 補間した結果を出力
    plt.plot(a, b, 'k', label=interpolation_method)
   
    # グラフのタイトル・目盛設定
    plt.title(interpolation_method)
    plt.grid(True)
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # 凡例の位置設定
    plt.legend(loc='lower right')
   
    # グラフ出力
    file_name = interpolation_method + '.jpg'
    plt.savefig(file_name)
   
    # 表示
    plt.show()
 
#==============================================================================
# 決定係数が最大を示す n 次式の最小二乗法近似関数
#==============================================================================
def n_th_order_least_squares_approximation_by_max_corrcoef(x, y, stop_n, point):
   
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    max_corrcoef = 0                                       # 最終的に採用した決定係数記憶用変数
    orders = []                                            # 次数の履歴記憶用変数
    corrcoefs = []                                         # 決定係数の履歴記憶用変数
    for n in range(1, stop_n):                             # 1 ~ stop_n まで繰り返す
        f = np.poly1d(np.polyfit(x, y, n))                 # 多項式近似
        Y = f(X)                                           # X に対応する y の値を計算する
        corrcoef = np.corrcoef(y, f(x))[0, 1]**2           # 決定係数を計算する
        orders.append(n)                                   # 次数の履歴記憶
        corrcoefs.append(corrcoef)                         # 決定係数の履歴記憶
        # 現状最大の決定係数 max_corrcoef よりさっき計算した corrcoef が大きい場合
        if max_corrcoef < corrcoef:
            max_corrcoef = corrcoef                        # max_corrcoef の更新
            fix_Y = Y                                      # 最大決定係数の時の Y を記憶する
            fix_n = n                                      # 最大決定係数を記憶する
    plt.figure(figsize=(8, 6))                             # グラフの初期化
    plt.plot(orders, corrcoefs, 'r.', label="coefficient of determination") # 決定係数の値の出力
    plt.plot(orders, corrcoefs, 'r')                       # 値のドットを線で結ぶ
    plt.legend(loc='lower right')                          # 凡例の位置設定
    plt.grid(True)                                         # 罫線表示
    plt.xlabel('order')                                    # 横軸タイトル
    plt.ylabel('coefficient of determination')             # 縦軸タイトル
    file_name = 'Determination coefficient transition to ' + str(fix_n) + '.jpg' # グラフ名設定
    plt.savefig(file_name)                                 # グラフ出力
    return X, fix_Y, fix_n                                 # 戻り値
 
#==============================================================================
# 決定係数が上昇しなくなった n 次式の最小二乗法近似関数
#==============================================================================
def n_th_order_least_squares_approximation_by_eps(x, y, stop_n, point):
   
    X = np.linspace(x[0], x[-1], num=point, endpoint=True) # x の値の準備
    max_corrcoef = 0                                       # 最終的に採用した決定係数記憶用変数
    orders = []                                            # 次数の履歴記憶用変数
    corrcoefs = []                                         # 決定係数の履歴記憶用変数
    eps = 0.00001                                          # 打切り誤差
    for n in range(1, stop_n):                             # 1 ~ stop_n まで繰り返す
        f = np.poly1d(np.polyfit(x, y, n))                 # 多項式近似
        Y = f(X)                                           # X に対応する y の値を計算する
        corrcoef = np.corrcoef(y, f(x))[0, 1]**2           # 決定係数を計算する
        orders.append(n)                                   # 次数の履歴記憶
        corrcoefs.append(corrcoef)                         # 決定係数の履歴記憶
        # 現状最大の決定係数 max_corrcoef とさっき計算した corrcoef の差の絶対値が打切り誤差より小さい場合
        if abs(max_corrcoef - corrcoef) < eps:
            break                                          # 多項式計算のループ終了
        # 現状最大の決定係数 max_corrcoef よりさっき計算した corrcoef が大きい場合
        if max_corrcoef < corrcoef:
            max_corrcoef = corrcoef                        # max_corrcoef の更新
            fix_Y = Y                                      # 最大決定係数の時の Y を記憶する
            fix_n = n                                      # 最大決定係数を記憶する
    plt.figure(figsize=(8, 6))                             # グラフの初期化
    plt.plot(orders, corrcoefs, 'r.', label="coefficient of determination") # 決定係数の値の出力
    plt.plot(orders, corrcoefs, 'r')                       # 値のドットを線で結ぶ
    plt.legend(loc='lower right')                          # 凡例の位置設定
    plt.grid(True)                                         # 罫線表示
    plt.xlabel('order')                                    # 横軸タイトル
    plt.ylabel('coefficient of determination')             # 縦軸タイトル
    file_name = 'Determination coefficient transition to ' + str(fix_n) + '.jpg' # グラフ名設定
    plt.savefig(file_name)                                 # グラフ出力
    return X, fix_Y, fix_n                                 # 戻り値
 
#==============================================================================
# テスト用データ作成
#==============================================================================
# 値の準備 =====================================================================
N = 2**7                           # サンプリング数
dt = 0.00025                       # サンプリング周期
ampl = np.array([ 1,  1,  2])      # 振幅
freq = np.array([60, 45, 100])     # 周波数
peri_axis = np.arange(0, N*dt, dt) # 周期軸
error = 0.15                       # 誤差係数
for i in range(len(freq)):
    if i==0:
        origin_y = ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis)                       # 入力波形
        y = ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis) + error * np.random.randn(N) # 誤差を与えた入力波形
    else:
        origin_y += ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis)                       # 入力波形
        y += ampl[i] * np.sin(2*np.pi*freq[i]*peri_axis) + error * np.random.randn(N) # 誤差を与えた入力波形
 
# 与えた誤差の計算 ==============================================================
err = (y) - (origin_y)
 
# グラフ出力設定 ================================================================
xmin = min(peri_axis) - (3*dt) # グラフの x 軸の最小値設定
xmax = max(peri_axis) + (3*dt) # グラフの x 軸の最大値設定
ymin = min(y) * 1.2            # グラフの y 軸の最小値設定
ymax = max(y) * 1.2            # グラフの y 軸の最大値設定
x_scale_interval = 1           # x 方向目盛間隔
y_scale_interval = 1           # y 方向目盛間隔
 
max_err = max(err)
min_err = min(err)
ymin_err = min_err * 1.5
ymax_err = max_err * 1.5
 
# グラフ出力 ====================================================================
wave_graph(peri_axis, origin_y, peri_axis, y, xmin, xmax, ymin, ymax)
err_graph(peri_axis, err, xmin, xmax, ymin_err, ymax_err)
 
# 決定係数が最大を示す次数による多項式補間 =========================================
X, Y, n = n_th_order_least_squares_approximation_by_max_corrcoef(peri_axis, y, 100, 1000)
graph(peri_axis, origin_y, X, Y, str(n) + 'th order by max coef')
 
# 決定係数が上昇しなくなった次数による多項式補間 ======================================
X, Y, n = n_th_order_least_squares_approximation_by_eps(peri_axis, y, 100, 1000)
graph(peri_axis, origin_y, X, Y, str(n) + 'th order by stop coef')



結果

曲線の近似には元々の波形と近似曲線がどれくらい一致しているかという指標があります。

この指標を"決定係数"と言います。
決定係数が示す範囲は 0 〜 1 です。
1 に近づく程、波形が一致していることを示します。

しかし、この指標は注意が必要で、波形データは点の集まりなので、その点における一致を見る指標だからです。
イメージがつかないと思うので、実際に結果を見てみます。

今回は色々と次数を 1 次から 2 次、 3 次と 100 次まで上げていきました。

次数が上がるにつれて決定係数も上昇していく傾向があります。

そこで、下記の①、②の近似曲線を作成して傾向を確認しました。

①決定係数が最大を示した次数の近似曲線
②決定係数の上昇が無くなった次数の近似曲線

図1 作成波形グラフ


図2 誤差グラフ


図3 ①の決定係数グラフ


図4 ①の近似曲線グラフ


図5 ②の決定係数グラフ


図6 ②の近似曲線グラフ


①の近似曲線は 99 次という大きな次数の近似曲線です。
図3 から次数を上げるにつれて決定係数も上昇している傾向が分かります。
決定係数が 1 付近に来るとほぼ一定になります。

図4 は 0.035 秒あたりで波形が乱れています。
決定係数は高く、一致しているかと思われましたが、点における一致を見る指標であるため、点と点の間の乱れまでは分かりません。

図5 を見るとほぼ一定になり、すぐに終了しています。

図6 は乱れがなく図4 と遜色ない程度に近似されています。

決定係数の値としては②の近似曲線よりも①の近似曲線の方が高いはずですが、②の近似曲線の方が曲線としては一致していると言えます。

決定係数は点がどのくらい一致するのかを示す指標なので、決定係数の高さだけでは最適な次数の近似曲線を決めることはできないことが分かりました。



コメント

今回、近似曲線を求めた時に下記の Warning がスローされました。


C:\WinPython-64bit-3.4.4.6Qt5\python-3.4.4.amd64\lib\site-packages\numpy\lib\polynomial.py:595: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

調べたら、高次の多項式近似の場合、近似の信頼性が低下するらしく、その時にスローされる Warning の様でした。

より高次の複雑な曲線の近似方法について調べていきたいです。



以上

Python で信号のノイズ除去 〜移動平均〜



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


この記事の目次



背景・目的


以前の記事で移動平均で波形データのスムージングが出来ることが分かりました。
Python で信号のノイズ除去 〜ローパスフィルター〜 - 解析エンジニアの自動化 blog

下のリンクのローパスフィルターと振幅処理でテストした波形データで試してみます。

Python で信号のノイズ除去 〜ローパスフィルター〜 - 解析エンジニアの自動化 blog

Python で信号のノイズ除去 〜振幅処理〜 - 解析エンジニアの自動化 blog



動作環境


Windows 7
・winpython 64bit 3.4.4



テスト用波形

テストに使用した波形データの仕様を以下に示す。

なお、テスト用波形データはリンクの記事の波形データ作成プログラムにノイズを与えられるように改造して作成しました。

Python で複雑な波形データを作る - 解析エンジニアの自動化 blog


正弦波数:1波
  サンプリング点数:1024点
  サンプリング周期:0.05秒
  正弦波式: A × sin( 2 × π × f × t )
 
正弦波式
 テスト用波形の正弦波の式を示す。
  0.869669428796 × sin( 2 × π × 0.0675797949804 × t )

Wave for Test



プログラム

ソースコード


###############################################################################
# 波形分析プログラム - Waveform analysis program -
###############################################################################
#==============================================================================
# ファイルインポート
#==============================================================================
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
 
#==============================================================================
# 波形デジタルデータの csv 入力関数
#==============================================================================
def ReadWaveCsv():
   
    #ファイル読み込み ############################################################
    wave = np.loadtxt('C:\\Users\\UserName\\Desktop\\test-wave.csv', delimiter=',')
   
    # wave データの列スライス #####################################################
    time = wave[0:, 0] # time[sec]
    ampl = wave[0:, 1] # amplitude
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(time, ampl)
    plt.grid(True)
    plt.title('Input Wave')
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Input_Wave.jpg'
    plt.savefig(file_name)
   
    return time, ampl, file_name
 
#==============================================================================
# 移動平均する関数
#==============================================================================
def moving_average(x, y, term):
   
    # カーネル作成 ##############################################################
    kernel = np.ones(term)/float(term)
   
    # 移動平均 ################################################################
    y2 = np.convolve(y, kernel, mode='valid')
    x2 = np.convolve(x, kernel, mode='valid')
   
    # グラフの作成 ##############################################################
    plt.figure(figsize=(10, 8))
    plt.plot(x,  y,  'r--', label='original wave')
    plt.plot(x2, y2, 'k-',  label='moving average term = ' + str(term))
    plt.grid(True)
    plt.title('Moving Average')
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
    plt.legend()
   
    # グラフの出力 ##############################################################
    file_name = 'smoothing(moving average term = ' + str(term) + ').jpg'
    plt.savefig(file_name)
   
    return x2, y2, file_name
 
#==============================================================================
# 実用例
#==============================================================================
# 波形の読込
time, ampl, input_wave_name = ReadWaveCsv()
 
# 移動平均
ret_time, ret_ampl, ma_name = moving_average(time, ampl, 20)
 
# グラフ表示
plt.figure(figsize=(10, 8))
plt.plot(time, ampl)
plt.plot(ret_time, ret_ampl)
plt.grid(True)
plt.title('Noise Cleared Wave and Original Wave')
plt.xlabel('time[sec]')
plt.ylabel('amplitude')
 
# グラフ出力
file_name = 'Noise_Cleared_Wave_and_Original_Wave.jpg'
plt.savefig(file_name)



結果

このプログラムから出力されるノイズ入り波形・移動平均波形・ノイズ除去波形とノイズ入り波形の重ね合わせをそれぞれ図1 〜 3 に示します。

図1 Input Wave


図2 Moving Average


図3 Noise Cleared Wave and Original Wave



コメント

単なる平均化なので、ノイズが小さくなってはいるもののノイズ除去波形自体がキレイではありません。

ローパスフィルターや振幅カットなどで処理した方が良さそうです。

移動平均は波形分析手法とはちょっと違う気がしました。



以上

ノイズを含む単純波形を移動平均でスムージングする



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


この記事の目次



背景・目的


前の記事で波形に移動平均でスムージングする方法を調べました。

移動平均の実装方法を調べただけで、移動平均によるスムージングの効果を確認していなかったので、テスト用の簡単な波形データを作り、 その波形に誤差を与えます。

そして、 元の簡単な波形に近づけてみたいと思います。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 移動平均によるスムージングプログラム
###############################################################################
# ファイルインポート ###############################################################
import numpy as np
import matplotlib.pyplot as plt
 
# テスト波形の作成 ###############################################################
x     = np.linspace(0, 10, 100)              # x 軸の値
y_org = np.cos(x)                            # y の値
np.random.seed(seed=12)
y     = np.cos(x) + np.random.randn(100)*0.1 # ノイズを含んだ y の値
 
for i in range(2, 10):
   
    # 移動平均 ################################################################
    term = i # 移動平均の個数
    kernel = np.ones(term)/term # カーネル
    x2 = np.convolve(x, kernel, mode='valid') # y の移動平均
    y2 = np.convolve(y, kernel, mode='valid') # x の移動平均
   
    # グラフの出力 ##############################################################
    plt.figure(figsize=(10, 8))
    plt.plot(x,  y_org, 'r',   label='original-sin')
    plt.plot(x,  y,     'k--', label='original-sin with noise')
    plt.plot(x2, y2,    'g-',  label='moving average term = ' + str(term))
    plt.legend()
   
    # グラフの出力 ##############################################################
    file_name = 'smoothing(moving average term = ' + str(term) + ').jpg'
    plt.savefig(file_name)



結果

このプログラムから出力される元波形・ノイズ入り波形・ノイズ除去波形の重ね合わせグラフを図1 に示します。

図1 グラフ term = 2


図2 グラフ term = 3


図3 グラフ term = 4


図4 グラフ term = 5


図5 グラフ term = 6


図6 グラフ term = 7


図7 グラフ term = 8
図8 グラフ term = 9



コメント

意外とキレイに除去されています。
複雑な波形に対して効果があるのか…

term による違いが見たかったので、与える誤差は seed 値を固定して、 term を変更して繰り返し実行しても同じ乱数が発生するようにしました。

term が大きくなるにつれて元波形に近づいてます。

しかし、 term を大きくしたり、何度も移動平均するとピーク値が小さくなるデメリットがありますね。
使い方が難しい‼︎



以上