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

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

Python で曲線近似(カーブフィッティング)する



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


この記事の目次



背景・目的


高速フーリエ変換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 次スプラインよりもさらに滑らかになっています。

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

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

図1 1 次スプライン補間


図2 2 次スプライン補間


図3 3 次スプライン補間


図4 秋間補間


図5 全補間重ね合わせ



コメント

波形の様なデータは 2 次、 3 次、秋間のスプライン補間を使えば良さそうな事が分かりました。



以上

Python で複雑な波形の分析入門



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


この記事の目次



背景・目的


Python で読み込んだ波形データを 高速フーリエ変換FFT )するプログラムを作成しました。

普通の FFT です。
特に何かがすごいということもありません。

読み込んだ波形データに対し、 FFT をして周波数、振幅を読み取って、 Word ファイルにまとめています。
さらに、 FFT と読み取った値の検証として、読み取った周波数、振幅の正弦波を重ね合わせて読み込んだ波形データと一致するかを行なっています。

あったら便利だし、使うとは思います。

勉強の為に分析結果を Microsoft Word で作成する機能をつけました。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 波形の分析に係る関数群
###############################################################################
 
#==============================================================================
# ファイルインポート
#==============================================================================
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import docx
from docx.enum.text import WD_ALIGN_PARAGRAPH
import math
 
#==============================================================================
# 波形デジタルデータの csv 入力関数
#==============================================================================
def ReadWaveCsv():
    #ファイル読み込み
    wave = np.loadtxt('C:\\Users\\UserName\\Desktop\\python-csv.csv', delimiter=',')
   
    time = wave[0:, 0] # time[sec]
    ampl = wave[0:, 1] # amplitude
   
    # グラフ表示
    plt.figure()
    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 fft(time, ampl, file_name):
    # サンプリング周期[sec]の計算
    dt = time[1] - time[0]
   
    # データ数カウント
    N = len(ampl)
   
    # 高速フーリエ変換( FFT )
    f = np.fft.fft(ampl)
   
    # FFT の複素数結果を絶対に変換
    absf = np.abs(f)
   
    # 振幅をもとの信号に揃える
    absf_amp = absf / N * 2       # 交流成分
    absf_amp[0] = absf_amp[0] / 2 # 直流成分
   
    # 周波数軸のデータ作成
    fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始, 終了, 分割数)
   
    # ピーク検出
    maximal_idx = signal.argrelmax(absf_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
    for i in range(len(maximal_idx[0])):
        print(str(fq[maximal_idx[0][i]]) + ' Hz, ' + str(absf_amp[maximal_idx[0][i]]))
   
    # グラフ表示
    plt.figure()
    plt.plot(fq[:int(N/2)+1], absf_amp[:int(N/2)+1])
    plt.scatter(fq[maximal_idx], absf_amp[maximal_idx], c='red', s=25)
    plt.grid(True)
    plt.title('Fast Fourier Transform')
    plt.xlabel('freqency[Hz]')
    plt.ylabel('amplitude')
   
    # グラフ出力
    file_name = 'Fast_Fourier_Transform.jpg'
    plt.savefig(file_name)
   
    return dt, N, fq, absf_amp, maximal_idx, file_name
 
#==============================================================================
# FFT で得られた正弦波式を合成して分析対象波形と比較検証する関数
#==============================================================================
def VerificationSynthesis(data_number, dt, maximal_idx, absf_amp, fq, time, ampl):
    t = np.arange(0, data_number*dt, dt)
   
    for i in range(len(maximal_idx[0])):
        if i==0:
            f  = absf_amp[maximal_idx[0][i]] * np.sin(2 * np.pi * fq[maximal_idx[0][i]] * t)
        else:
            f += absf_amp[maximal_idx[0][i]] * np.sin(2 * np.pi * fq[maximal_idx[0][i]] * t)
   
    # グラフ表示
    plt.figure()
    plt.plot(time, ampl)
    plt.plot(t, f)
    plt.grid(True)
    plt.title('Verification Synthesis Wave on Input Wave')
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力
    file_name = 'Verification.jpg'
    plt.savefig(file_name)
   
    return f, file_name
 
#==============================================================================
# データ間の類似度を計算する関数(ユークリッド距離)
#==============================================================================
def EuclidDistance(signal1, signal2):
    dis = signal2 - signal1
    disdis = dis*dis
    euclid = math.sqrt(np.sum(disdis))
    print(euclid)
    return euclid
 
#==============================================================================
# 試験用波形を作成する関数 Wave の実行結果報告書作成関数
#==============================================================================
def DocumentReport(time, ampl, dt, N, input_wave_name, fft_name, maximal_idx, fq, absf_amp, verification_name, euclid):
   
    # add_heading の第2引数
    # 0 : 表題
    # 1 : 見出し1
    # 2 : 見出し2
    # 3 : 見出し3
    # 4 : 見出し4
   
    # Microsoft Word オブジェクト生成
    doc = docx.Document()
   
    # 見出し:表題の追加
    doc.add_heading('波形の分析結果', 0)
   
    # 見出し:見出し1
    doc.add_heading('波形の分析方法', 1)
    doc.add_paragraph(' 高速フーリエ変換を実施し、波形を構成する単純正弦波を求める。\n')
   
    # 見出し:見出し1
    doc.add_heading('分析対象波形', 1)
    doc.add_paragraph(' 分析した波形を図1に示す。また、波形のデジタルデータを添付資料-1に示す。')
    doc.add_picture(input_wave_name) # 作成した波形グラフの挿入
    p = doc.add_paragraph('図1 分析した波形\n')
    p.alignment = WD_ALIGN_PARAGRAPH.CENTER
    
    # 見出し:見出し1
    doc.add_heading('分析結果', 1)
    doc.add_paragraph(' 高速フーリエ変換結果を図2に示し、波形仕様及び波形を構成する正弦波式をまとめる。')
    doc.add_picture(fft_name) # 作成した波形グラフの挿入
    p = doc.add_paragraph('図2 高速フーリエ変換結果\n')
    p.alignment = WD_ALIGN_PARAGRAPH.CENTER
   
    # 見出し:見出し2
    doc.add_heading('波形仕様', 2)
    doc.add_paragraph(' 波形仕様を列挙する。')
    doc.add_paragraph('  合成正弦波数:' + str(len(maximal_idx[0])) + '波')
    doc.add_paragraph('  サンプリング点数:' + str(N) + '点')
    doc.add_paragraph('  サンプリング周期:' + str(dt) + '秒\n')
 
    # 見出し:見出し2
    doc.add_heading('正弦波式', 2)
    doc.add_paragraph(' 高速フーリエ変換結果から正弦波式を列挙する。')
    for i in range(len(maximal_idx[0])):
        doc.add_paragraph('  ' + str(absf_amp[maximal_idx[0][i]]) + ' × sin( 2 × π × ' + str(fq[maximal_idx[0][i]]) + ' t )')
    doc.add_paragraph('')
   
    # 見出し:見出し1
    doc.add_heading('検証', 1)
    doc.add_paragraph(' 高速フーリエ変換で得られた正弦波式を分析対象波形と同仕様で合成し、分析対象波形と比較する。')
    doc.add_picture(verification_name)
    p = doc.add_paragraph('図3 高速フーリエ変換によって得た正弦波式の合成波形\n')
    p.alignment = WD_ALIGN_PARAGRAPH.CENTER
    if euclid < 30:
        doc.add_paragraph(' 図3より、分析対象波形と高速フーリエ変換で得られた正弦波式の重ね合わせ波形はその大部分が再現出来ている。\n')
    elif euclid < 45:
        doc.add_paragraph(' 図3より、分析対象波形と高速フーリエ変換で得られた正弦波式の重ね合わせ波形は類似傾向にはあるが誤差は大きい。\n')
    else:
        doc.add_paragraph(' 図3より、誤差が大きく高速フーリエ変換で得られた正弦波式の重ね合わせ波形では分析対象波形の再現は出来なかった。\n')
   
    # 見出し:表題の追加
    doc.add_heading('添付資料-1', 0)
    doc.add_paragraph('分析対象波形デジタルデータ')
    doc.add_paragraph('time[sec],amplitude')
    for i in range(len(ampl)):
        doc.add_paragraph(str(time[i]) + ',' + str(ampl[i]))
   
    # 保存
    doc.save('C:\\Users\\UserName\\Desktop\\Analysis_Report.docx')
 
#==============================================================================
# 実用例
#==============================================================================
# 波形の読込
time, ampl, input_wave_name = ReadWaveCsv()
 
# 高速フーリエ変換
dt, N, fq, absf_amp, maximal_idx, fft_name = fft(time, ampl, input_wave_name)
 
# 高速フーリエ変換によって得られた正弦波を重ね合わせ波形と分析対象波形の比較
f, verification_name = VerificationSynthesis(N, dt, maximal_idx, absf_amp, fq, time, ampl)
 
# データ間の類似度計算
euclid = EuclidDistance(ampl, f)
 
# Misrosoft Word での波形分析レポート作成
DocumentReport(time, ampl, dt, N, input_wave_name, fft_name, maximal_idx, fq, absf_amp, verification_name, euclid)



結果

FFT から得られた周波数、振幅の正弦波を重ね合わせて読み込んだ波形データと一致するかを試したのですが、完全には一致しませんでした。

あと、波形データ間の類似度計算はユークリッド距離にしています。

判定も "DocumentReport" 関数の中でかなりテキトーにやっているので、信頼できない類似度です。

何か良い類似度計算方法を思いついたら、また記事にしたいと思います。

図1 読み込んだ波形データ


図2 FFT 結果


図3 入力波形と検証波形の比較


図4 Word ファイル(冒頭)



コメント

FFT から得られた周波数、振幅の波形の重ね合わせが読み込んだ波形と一致しなくて残念でした。

FFT の精度向上について調べていきたいです。



以上

Python で複雑な波形データを作る



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


この記事の目次



背景・目的


Python で分析をするにあたり、今のところ波形データの分析が取り組みやすく、簡単だなぁ…と思っているところです。

分析対象の波形作成を簡単かつスピーディにするためのプログラムを作りました。

Microsoft Word で作成した波形の仕様書もついでに出力するプログラムにしました。

Microsoft Word の操作は下の記事でまとめています。
Python で docx を使って Microsoft Word を操作する - 解析エンジニアの自動化 blog



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 試験用波形の作成に係る関数群
###############################################################################
 
#==============================================================================
# ファイルインポート
#==============================================================================
import numpy as np
import matplotlib.pyplot as plt
import docx
from docx.enum.text import WD_ALIGN_PARAGRAPH
 
#==============================================================================
# 試験用波形を作成する関数
#==============================================================================
def Wave(wave_number, data_number, dt):
   
    ampl        = np.random.rand(wave_number) # 0.0 以上 1.0 未満の乱数
    freq        = np.random.rand(wave_number) # 0.0 以上 1.0 未満の乱数
    t           = np.arange(0, data_number*dt, dt)
   
    for i in range(0, wave_number):
        if i==0:
            f  = ampl[i] * np.sin(2 * np.pi * freq[i] * t)
        else:
            f += ampl[i] * np.sin(2 * np.pi * freq[i] * t)
   
    # グラフ表示
    plt.plot(t, f)
    plt.grid(True)
    plt.title('Wave')
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力
    file_name = 'wave.jpg'
    plt.savefig(file_name)
   
    return ampl, freq, t, f, file_name
 
#==============================================================================
# 試験用波形を作成する関数 Wave の実行結果報告書作成関数
#==============================================================================
def DocumentReport(wave_number, data_number, dt, ampl, freq, t, f, img_name):
   
    # add_heading の第2引数
    # 0 : 表題
    # 1 : 見出し1
    # 2 : 見出し2
    # 3 : 見出し3
    # 4 : 見出し4
   
    # Microsoft Word オブジェクト生成
    doc = docx.Document()
   
    # 見出し:表題の追加
    doc.add_heading('試験用波形の作成', 0)
   
    # 見出し:見出し1
    doc.add_heading('試験用波形の作成方法', 1)
    doc.add_paragraph(' 試験用波形は正弦波の重ね合わせによって作成する。\n')
   
    # 見出し:見出し1
    doc.add_heading('試験用波形の作成仕様', 1)
    doc.add_paragraph(' 試験用波形作成仕様を列挙する。\n')
    doc.add_paragraph('  合成正弦波数:' + str(wave_number) + '波')
    doc.add_paragraph('  サンプリング点数:' + str(data_number) + '点')
    doc.add_paragraph('  サンプリング周期:' + str(dt) + '秒')
    doc.add_paragraph('  正弦波式: A × sin( 2 × π × f × t )\n')
   
    # 見出し:見出し2
    doc.add_heading('正弦波式', 2)
    doc.add_paragraph(' 試験用波形の正弦波の重ね合わせに使用した正弦波式を列挙する。')
    for i in range(len(ampl)):
        doc.add_paragraph('  ' + str(ampl[i]) + ' × sin( 2 × π × ' + str(freq[i]) + ' × t )')
    doc.add_paragraph('')
   
    # 見出し:見出し1
    doc.add_heading('作成結果', 1)
    doc.add_paragraph(' 作成した試験用波形を図1に示し、試験用波形のデジタルデータをまとめる。')
    doc.add_picture(img_name) # 作成した波形グラフの挿入
    p = doc.add_paragraph('図1 試験用波形\n')
    p.alignment = WD_ALIGN_PARAGRAPH.CENTER
   
    doc.add_page_break() # 改ページ
   
    # 見出し:見出し1
    doc.add_heading('デジタルデータ', 0)
    doc.add_paragraph('time[sec],amplitude')
    for i in range(len(f)):
        doc.add_paragraph(str(t[i]) + ',' + str(f[i]))
   
    # 保存
    doc.save('C:\\Users\\UserName\\Desktop\\python-docx.docx')
 
#==============================================================================
# 試験用波形を作成する関数 Wave で作成した試験用波形デジタルデータの csv 出力関数
#==============================================================================
def WriteWaveCsv(t, f):
    w = np.stack([t, f])
    np.savetxt('C:\\Users\\UseName\\Desktop\\python-csv.csv', w.T, delimiter=',')
 
#==============================================================================
# 実用例
#==============================================================================
wave_number = 5     # 合成する波形の数
data_number = 2**12 # データ数
dt = 0.01           # サンプリング周期 [sec]
 
# 波形作成
ampl, freq, t, f, img_name = Wave(wave_number, data_number, dt)
 
# Misrosoft Word での波形作成レポート作成
DocumentReport(wave_number, data_number, dt, ampl, freq, t, f, img_name)
 
# 波形データの csv 出力
WriteWaveCsv(t, f)



結果

このプログラムから出力される波形のグラフ、 Word ファイル、 csv ファイルをそれぞれ図1 〜 図3 に示します。

図2 と図3 の Word ファイルと csv ファイルは内容が1枚の画像には収まらないので、冒頭の部分のみの画像にしました。

図1 作成した波形


図2 作成した Word ファイル(冒頭)


図3 作成した csv ファイル(冒頭)



コメント

こういったプログラムを量産していくと仕事の自動化に効果的な気がしてます。

さらに、 GUI 化や exe 化などで汎用的に出来れば効果はもっと大きいものとなると思います。

GUI の作り方と exe 化については下の記事でまとめています。

GUI の作り方】
Python で ファイル選択ダイアログを使う - 解析エンジニアの自動化 blog

【 exe 化】
Python を exe ファイルに変換する - 解析エンジニアの自動化 blog



以上

Python で docx を使って Microsoft Word を操作する



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


この記事の目次



背景・目的


Python に仕事をさせて最後に待ち受けているのが報告書の作成です。

自動化している事を悟られない様に Microsoft Word で報告書を書きました。

docx というモジュールの入門みたいな記事です。



動作環境


Windows 7
・winpython 64bit 3.4.4



インストール

pip コマンドからインストールしました。

pip install python-docx



プログラム

ソースコード


import docx
from docx.enum.text import WD_ALIGN_PARAGRAPH

# ワードドキュメント作成
doc = docx.Document()

# 段落に見出し追加
doc.add_heading('報告書', 0)

# 段落に文章追加
doc.add_heading('はじめに', 1)
doc.add_paragraph('はじめに〜')

doc.add_heading('対象', 1)
doc.add_paragraph('対象〜')
p = doc.add_paragraph('図1 対象')
p.alignment = WD_ALIGN_PARAGRAPH.CENTER

doc.add_heading('方法', 1)
doc.add_paragraph('方法〜')

doc.add_heading('結果', 1)
doc.add_paragraph('結果〜')
p = doc.add_paragraph('図2 結果')
p.alignment = WD_ALIGN_PARAGRAPH.CENTER

# ファイル名を指定して保存
doc.save('C:\\Users\\UserName\\Desktop\\test-docx.docx')



結果

デスクトップに "test-docx.docx" というファイルが出来たと思います。

図1 に作成した Word ファイルを示します。

図1 作成した Word ファイル



コメント

Python から Microsoft-Word の操作は簡単です!

目次やヘッダー・フッターなど、 Microsoft-Word で使う機能は言うほど多くないので、もう少し調べれば普段の仕事で行う Microsoft-Word のタスクは自動化出来そうです。

以外と、インターネットに情報が少なくて驚きました。

ソースコードの分割が出来れば Word ファイルを操作する Python 関数を量産したソースコードをひとまとめに出来ます。

ソースコードの分割は下記の記事にサンプルをまとめています。
Python のソースコードを分割する - 解析エンジニアの自動化 blog



以上

Python のソースコードを分割する



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


この記事の目次



背景・目的


Python でコードを書いているときに、クラス毎や関数毎にファイルが分割出来たら良いなと思いました。
クラス毎や関数毎に分担する事も出来そうですし…

ファイルが分割出来ると1度作ったソースコードの流用がとても簡単になり、ソースコードを資産的に管理する事が出来る様になります。

なので、ファイル分割の方法をまとめます。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコードは以下の2つを用意します。
関数とクラスのファイル:split_files_child.py
関数とクラスの実行ファイル:split_files.py

ソースコード 〜関数とクラスのファイル〜


###############################################################################
# サンプル関数・クラスファイル
###############################################################################
 
# 関数 1
def sample_function1():
    print('split_files_child.py の sample_function1 が実行されました。')
 
# 関数 2
def sample_function2():
    print('split_files_child.py の sample_function2 が実行されました。')
 
# クラス
class samplclass():
   
    # コンストラクタ
    def __init__(self, arg1):
        self.value = arg1
   
    # メソッド
    def method1(self):
        self.value = '\n' + self.value + '\nsampleclass の method1 によって追加された文章!!'



ソースコード 〜別ファイルの関数とクラスの実行ファイル〜


###############################################################################
# 分割されているソースコードファイルの実行プログラム
###############################################################################
 
# ファイルインポート
import split_files_child as sfc # .py は不要
 
# split_files_child.py の sample_function1 を実行
sfc.sample_function1()
 
# split_files_child.py の sample_function1 を実行
sfc.sample_function2()
 
# split_files_child.py の sampleclass のインスタンス生成
# 'split_files.py' で初期化
myclass = sfc.samplclass('split_files.py')
 
# value の表示
print(myclass.value)
 
# split_files_child.py の sampleclass のメソッド実行
myclass.method1()
 
# value の表示
print(myclass.value)



結果

かなり簡単に分割出来ました。

関数とクラスのファイルを実行ファイルと同じディレクトリに置いて、 import しただけです。


In [1]: runfile('C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3/split_files.py', wdir='C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3')
Reloaded modules: split_files_child
split_files_child.py の sample_function1 が実行されました。
split_files_child.py の sample_function2 が実行されました。
split_files.py
 
split_files.py
sampleclass の method1 によって追加された文章!!



コメント

インターネットを調べているとクラスだけ書かれているサイトを見かけたりします。

こうやって使えば良かったんですね。

これで、コーディング時間が削減できたら良いですね。

あまり視覚的では無いので読みづらい記事になってしまった…



以上

Python で逆高速フーリエ変換(IFFT)をする



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


この記事の目次



背景・目的


フーリエ変換をやったので、逆フーリエ変換をやります。



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 逆高速フーリエ変換( IFFT )を計算するプログラム
###############################################################################
# インポート
import matplotlib.pyplot as plt
import numpy as np
 
# 簡単な信号の作成
N     = 2**6 # サンプル数
dt    = 0.025 # サンプリング周期( sec )
freq1 = 10 # 周波数( Hz )
ampl1 = 1 # 振幅
freq2 = 15 # 周波数( Hz )
ampl2 = 1 # 振幅
print("■■■ 入力条件 ■■■")
print("サンプル数 : " + str(N))
print("サンプリング周期 : " + str(dt) + ' sec')
print("サンプリング周波数 : " + str(1/dt) + ' Hz')
print("入力波1")
print("  周波数 : " + str(freq1) + ' Hz')
print("  振 幅 : " + str(ampl1))
print("入力波2")
print("  周波数 : " + str(freq2) + ' Hz')
print("  振 幅 : " + str(ampl2))
 
# 時間軸
t = np.arange(0, N*dt, dt)
print("サンプリング時間 : " + str(N*dt-dt) + ' sec')
print("サンプリング時刻 : " + str(t))
 
# {周波数 freq1, 振幅 amp1 の正弦入力波} + {周波数 freq2, 振幅 amp2 の正弦入力波}
f = ampl1*np.sin(2*np.pi*freq1*t) + ampl2*np.sin(2*np.pi*freq2*t)
# 周波数 freq1, 振幅 amp1 の正弦入力波
#f = amp1 * np.sin(2*np.pi*freq1*t)
 
# 高速フーリエ変換( FFT )
F = np.fft.fft(f)
 
# FFT の複素数結果を絶対に変換
absf = np.abs(F)
 
# 振幅をもとの信号に揃える
absf_amp = absf / N * 2
absf_amp[0] = absf_amp[0] / 2
 
# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸 linspace(開始, 終了, 分割数)
 
idx = np.argmax(f)
print("\n■■■ 入力信号特性 ■■■")
print("入力信号の最大idx : " + str(idx))
print("入力信号の最大振幅 : " + str(f[idx]))
print("入力信号の最大時刻 : " + str(t[idx]))
 
idx = np.array(absf_amp[:int(N/2)+1]) # コピー
idx = idx.argsort()[::-1] # 降順に並べ替えた時のインデックスを取得する
print("\n■■■ フーリエ変換信号特性 ■■■")
print("フーリエ変換信号の最大idx : " + str(idx[0]))
print("フーリエ変換信号の最大振幅 : " + str(absf_amp[idx[0]]))
print("フーリエ変換信号の最大周波数 : " + str(fq[idx[0]]))
 
print("\nフーリエ変換信号の次点idx : " + str(idx[1]))
print("フーリエ変換信号の次点振幅 : " + str(absf_amp[idx[1]]))
print("フーリエ変換信号の次点周波数 : " + str(fq[idx[1]]))
 
# 逆フーリエ変換( IFFT )
F_ifft = np.fft.ifft(F)
 
# 実数部の取得
F_ifft_real = F_ifft.real
 
# スライス
F_ifft_real[:10]
 
idx = np.argmax(F_ifft_real)
print("\n■■■ 逆フーリエ変換信号特性 ■■■")
print("逆フーリエ変換信号の最大idx : " + str(idx))
print("逆フーリエ変換信号の最大振幅 : " + str(F_ifft_real[idx]))
print("逆フーリエ変換信号の最大時刻 : " + str(t[idx]))
 
# グラフ表示
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 12))
 
# 信号のグラフ(時間軸)
axes[0, 0].plot(t, f)
axes[0, 0].set_title('Input Wave')
axes[0, 0].set_xlabel('time[sec]')
axes[0, 0].set_ylabel('amplitude')
axes[0, 0].grid(True)
 
# FFTのグラフ(周波数軸)
axes[0, 1].plot(fq[:int(N/2)+1], absf_amp[:int(N/2)+1])
axes[0, 1].set_title('Fast Fourier Transform')
axes[0, 1].set_xlabel('freqency[Hz]')
axes[0, 1].set_ylabel('amplitude')
axes[0, 1].grid(True)
 
# IFFTのグラフ(時間軸)
axes[1, 0].plot(t, F_ifft_real, c="g")
axes[1, 0].set_title('Inverse Fast Fourier Transform')
axes[1, 0].set_xlabel('time[sec]')
axes[1, 0].set_ylabel('amplitude')
axes[1, 0].grid(True)
 
# Input Wave と IFFT Wave の重ね合わせ
axes[1, 1].plot(t, f, c="g")
axes[1, 1].plot(t, F_ifft_real, c="b")
axes[1, 1].set_title('Input Wave and IFFT')
axes[1, 1].set_xlabel('time[sec]')
axes[1, 1].set_ylabel('amplitude')
axes[1, 1].grid(True)
 
# グラフの出力
file_dir  = 'C:\\Users\\UserName\\Desktop\\'
file_name = 'ifft'
fig.savefig(file_dir + file_name + '0.0.jpg', bbox_unches="tight")



結果

図1 の 4 つのグラフの説明をします。

Input Wave は入力波のグラフです。

Fast Fourier Transform は Input Wave を高速フーリエ変換したグラフです。

Inverse Fast Fourier Transform は Fast Fourier Transform を逆高速フーリエ変換したグラフです。

Input Wave and IFFT は Input Wave と Inverse Fast Fourier Transform のグラフを重ね合わせたグラフです。

図1 変換結果


Input Wave and IFFT の図から Fast Fourier Transform の波形を逆高速フーリエ変換して求めた波形が入力した Input Wave とピッタリ重なっています。

そのため、正しく逆高速フーリエ変換が出来た事がわかります。

Python の 出力結果を示します。

In [1]: runfile('C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3/fft-0.0.py', wdir='C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3')
■■■ 入力条件 ■■■
サンプル数 : 64
サンプリング周期 : 0.025 sec
サンプリング周波数 : 40.0 Hz
入力波1
  周波数 : 10 Hz
  振 幅 : 1
入力波2
  周波数 : 15 Hz
  振 幅 : 1
サンプリング時間 : 1.5750000000000002 sec
サンプリング時刻 : [ 0. 0.025 0.05 ..., 1.525 1.55 1.575]
 
■■■ 入力信号特性 ■■■
入力信号の最大idx : 57
入力信号の最大振幅 : 1.70710678119
入力信号の最大時刻 : 1.425
 
■■■ フーリエ変換信号特性 ■■■
フーリエ変換信号の最大idx : 24
フーリエ変換信号の最大振幅 : 1.0
フーリエ変換信号の最大周波数 : 15.2380952381
 
フーリエ変換信号の次点idx : 16
フーリエ変換信号の次点振幅 : 1.0
フーリエ変換信号の次点周波数 : 10.1587301587
 
■■■ 逆フーリエ変換信号特性 ■■■
逆フーリエ変換信号の最大idx : 57
逆フーリエ変換信号の最大振幅 : 1.70710678119
逆フーリエ変換信号の最大時刻 : 1.425



コメント

Input Wave を FFT した後、ノイズを除去して、 IFFT する事があります。

Input Wave は横軸が時間のグラフで FFT したグラフは横軸が周波数です。

横軸が周波数になるとノイズなどが見やすくなることが多いので、波形をキレイにするときは、『入力波形 → FFT → ノイズ除去 → IFFT → キレイな波形』という手法は王道です。



以上

Python で高速フーリエ変換(FFT)をする



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


この記事の目次



背景・目的


とにかく何か分析めいたものがやりたくて、特に理由はないけれどフーリエ変換をやりました。

波形の分析の初歩なので、まとめておいて損はないと思いまして…



動作環境


Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 高速フーリエ変換( FFT )を計算するプログラム
###############################################################################
# インポート
import matplotlib.pyplot as plt
import numpy as np
 
# 簡単な入力波の作成
N     = 2**6 # サンプル数
dt    = 0.025 # サンプリング周期( sec )
freq1 = 10 # 周波数( Hz )
ampl1 = 1 # 振幅
freq2 = 15 # 周波数( Hz )
ampl2 = 1 # 振幅
print("■■■ 入力条件 ■■■")
print("サンプル数 : " + str(N))
print("サンプリング周期 : " + str(dt) + ' sec')
print("サンプリング周波数 : " + str(1/dt) + ' Hz')
print("入力波1")
print("  周波数 : " + str(freq1) + ' Hz')
print("  振 幅 : " + str(ampl1))
print("入力波2")
print("  周波数 : " + str(freq2) + ' Hz')
print("  振 幅 : " + str(ampl2))
 
# 時間軸
t = np.arange(0, N*dt, dt)
print("サンプリング時間 : " + str(N*dt-dt) + ' sec')
print("サンプリング時刻 : " + str(t))
 
# {周波数 freq1, 振幅 amp1 の正弦入力波} + {周波数 freq2, 振幅 amp2 の正弦入力波}
f = ampl1*np.sin(2*np.pi*freq1*t) + ampl2*np.sin(2*np.pi*freq2*t)
# 周波数 freq1, 振幅 amp1 の正弦入力波
#f = amp1 * np.sin(2*np.pi*freq1*t)
 
# 高速フーリエ変換( FFT )
F = np.fft.fft(f)
 
# FFT の複素数結果を絶対に変換
absf = np.abs(F)
 
# 振幅をもとの信号に揃える
absf_amp = absf / N * 2
absf_amp[0] = absf_amp[0] / 2
 
# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N)
 
idx = np.argmax(f)
print("\n■■■ 入力信号特性 ■■■")
print("入力信号の最大idx : " + str(idx))
print("入力信号の最大振幅 : " + str(f[idx]))
print("入力信号の最大時刻 : " + str(t[idx]))
 
idx = np.array(absf_amp[:int(N/2)+1]) # コピー
idx = idx.argsort()[::-1] # 降順に並べ替えた時のインデックスを取得する
print("\n■■■ フーリエ変換信号特性 ■■■")
print("フーリエ変換信号の最大idx : " + str(idx[0]))
print("フーリエ変換信号の最大振幅 : " + str(absf_amp[idx[0]]))
print("フーリエ変換信号の最大周波数 : " + str(fq[idx[0]]))
 
print("\nフーリエ変換信号の次点idx : " + str(idx[1]))
print("フーリエ変換信号の次点振幅 : " + str(absf_amp[idx[1]]))
print("フーリエ変換信号の次点周波数 : " + str(fq[idx[1]]))
 
# グラフ表示
fig, (axL, axR) = plt.subplots(ncols=2, figsize=(10,4))
 
# 信号のグラフ(時間軸)
axL.plot(t, f)
axL.set_title('Input Wave')
axL.set_xlabel('time[sec]')
axL.set_ylabel('amplitude')
axL.grid(True)
 
# FFTのグラフ(周波数軸)
axR.plot(fq[:int(N/2)+1], absf_amp[:int(N/2)+1])
axR.set_title('Fast ')
axR.set_xlabel('freqency[Hz]')
axR.set_ylabel('amplitude')
axR.grid(True)
 
# グラフの出力
file_dir  = 'C:\\Users\\UserName\\Desktop\\'
file_name = 'fft'
fig.savefig(file_dir + file_name + '0.0.jpg', bbox_unches="tight")



結果

numpy には高速フーリエ変換のメソッドがあります。
numpy は便利ですね。

図1 の Input Wave は分析した波形で、周波数 10Hz, 振幅 1 の正弦波と周波数 15Hz, 振幅 1 の正弦波の合成波でした。

図1 の Fast Fourier Transform は分析結果で、周波数が 10Hz と 15Hz でところで波形がとんがっています。
さらに、振幅がそれぞれ 1 を示しています。

見事に入力波が周波数 10Hz, 振幅 1 の正弦波と周波数 15Hz, 振幅 1 の正弦波の合成波である事が分かりました。

図1 変換結果



In [1]: runfile('C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3/fft0.0.py', wdir='C:/WinPython-64bit-3.4.4.6Qt5/settings/.spyder-py3')
■■■ 入力条件 ■■■
サンプル数 : 64
サンプリング周期 : 0.025 sec
サンプリング周波数 : 40.0 Hz
入力波1
  周波数 : 10 Hz
  振 幅 : 1
入力波2
  周波数 : 15 Hz
  振 幅 : 1
サンプリング時間 : 1.5750000000000002 sec
サンプリング時刻 : [ 0. 0.025 0.05 ..., 1.525 1.55 1.575]
 
■■■ 入力信号特性 ■■■
入力信号の最大idx : 57
入力信号の最大振幅 : 1.70710678119
入力信号の最大時刻 : 1.425
 
■■■ フーリエ変換信号特性 ■■■
フーリエ変換信号の最大idx : 24
フーリエ変換信号の最大振幅 : 1.0
フーリエ変換信号の最大周波数 : 15.2380952381
 
フーリエ変換信号の次点idx : 16
フーリエ変換信号の次点振幅 : 1.0
フーリエ変換信号の次点周波数 : 10.1587301587



コメント

入力波形がどんな単純波形の合成で出来ているのかが分かる様になりました。

自動車の振動の測定結果や観測した地震波がどんな単純波形の合成で出来ているのかがわかる様になるため、対策すれば良い単純な波形がわかる様になります。

ただし、振動の測定結果や観測地震波などは波形が複雑なため、構成される単純波形がものすごい数になる事でしょう…



以上