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

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

移動平均で波形のスムージングが出来る⁉︎ 〜 Python 〜



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


この記事の目次



背景・目的


フーリエ変換して周波数特性・振幅特性を利用してノイズ除去するとデータが振動したり、キレイにノイズ除去出来ない事もあります。

私もノイズを除去して波形が振動したので、何か良い方法はないかと調べていると移動平均でスムージングが出来るとの事でした。
気がつきませんでした…

Numpy の convolve メソッドで簡単に移動平均出来るとのことなので試してみたところ、かなり簡単に出来ました。

ただ、移動平均すると波形の両端が平均値を計算する関係上少なくなるので、グラフを描くときにつまずきました。

気づかなかったし、つまずいたのでまとめておきます。



動作環境


・Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# convolve メソッドを使用した移動平均プログラム
###############################################################################
# ファイルインポート ###############################################################
import numpy as np
import matplotlib.pyplot as plt
 
# テスト波形の作成 ###############################################################
#             1  2  3  4  5  6  7  8  9    要素数は 9 個
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
y = np.array([1, 3, 9, 3, 5, 5, 6, 1, 4])
 
for i in range(2, 9):
    # 移動平均 ################################################################
    term = i # 移動平均の個数
    b = np.ones(term)/term
    y2 = np.convolve(y, b, mode='valid') # x の移動平均
    x2 = np.convolve(x, b, mode='valid') # y の移動平均
   
    # グラフの作成 ##############################################################
    plt.figure(figsize=(8, 6))
    plt.plot(x,  y,  'k-',  label='original-sin with noise')
    plt.plot(x2, y2, 'r--', label='moving average term = ' + str(term))
    plt.legend()
    # グラフ出力
    file_name = 'moving average term = ' + str(term) + '.jpg'
    plt.savefig(file_name)



結果

このプログラムから出力される移動平均のグラフを図1 〜 図8 に示します。

図1 グラフ term = 2


図2 グラフ term = 3


図3 グラフ term = 4


図4 グラフ term = 5


図5 グラフ term = 6


図6 グラフ term = 7


図7 グラフ term = 8



コメント

グラフを描く上で x の値に手をこまねいた瞬間がありました。

結局、少し考えて x も移動平均してしまえば良いということに気づいて良かったです。

このスムージングにどの程度効果があるのか気になるところです。



以上

Python で信号のノイズ除去 〜振幅処理〜



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


この記事の目次



背景・目的


以前の記事でローパスフィルターをやりました。
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 fft(time, ampl):
   
    # サンプリング周期[sec]の計算 #################################################
    sampling_cycle = time[1] - time[0]
   
    # データ数カウント ############################################################
    N = len(ampl)
   
    # 高速フーリエ変換（ FFT ） ####################################################
    fft_ampl = np.fft.fft(ampl)
   
    # FFT の複素数結果を絶対に変換 ###############################################
    abs_fft_amp = np.abs(fft_ampl)
   
    # 振幅をもとの信号に揃える #####################################################
    abs_fft_amp    = abs_fft_amp    / N * 2   # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2       # 直流成分
   
    # 周波数軸のデータ作成 #######################################################
    frequency = np.linspace(0, 1.0/sampling_cycle, N) # 周波数軸　linspace(開始, 終了, 分割数)
   
    # ピーク検出 ################################################################
    maximal_idx = signal.argrelmax(abs_fft_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(frequency[:int(N/2)+1], abs_fft_amp[:int(N/2)+1])
    plt.scatter(frequency[maximal_idx], abs_fft_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 frequency, fft_ampl, file_name
 
#==============================================================================
# オペレーション・アンプリファイヤー 通称:オペアンプ
# cut_off 以下の amplitude をゼロにする
#==============================================================================
def operation_amplifier(fft_time, fft_amp, cut_off):
   
    # データ数カウント ############################################################
    N = len(fft_amp)
   
    # cut_off 処理用に実波形データに変換 ##########################################
   
    # FFT の複素数結果を絶対に変換
    abs_fft_amp = np.abs(fft_amp)
   
    # 振幅をもとの信号に揃える
    abs_fft_amp    = abs_fft_amp    / N * 2 # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2     # 直流成分
   
    # cut_off 以下の amplitude をゼロにする #######################################
    fft_amp[(abs_fft_amp <= cut_off)] = 0
   
    # グラフ用に実波形データに変換 #################################################
   
    # FFT の複素数結果を絶対に変換
    abs_fft_amp = np.abs(fft_amp)
   
    # 振幅をもとの信号に揃える
    abs_fft_amp    = abs_fft_amp    / N * 2 # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2     # 直流成分
   
    # ピーク検出 ################################################################
    #maximal_idx = signal.argrelmax(abs_fft_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
    maximal_idx = signal.argrelmax(abs_fft_amp, order=1) # 極大値インデックスの取得
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    #plt.plot(fft_time[:int(N/2)+1], abs_fft_amp[:int(N/2)+1])
    #plt.scatter(fft_time[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.plot(fft_time, abs_fft_amp)
    plt.scatter(fft_time[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.grid(True)
    plt.title('Operation Amplifier')
    plt.xlabel('freqency[Hz]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Operation_Amplifier.jpg'
    plt.savefig(file_name)
   
    return fft_time, fft_amp, file_name
 
#==============================================================================
# 波形を逆高速フーリエ変換する関数
#==============================================================================
def ifft(fft_time, fft_ampl):
   
    # データ数カウント ############################################################
    N = len(fft_ampl)
   
    # 逆高速フーリエ変換
    ampl = np.fft.ifft(fft_ampl)
   
    # 実数部分のみ
    ampl_real = ampl.real
   
    # 時間軸のデータ作成 ########################################################
    frequency = fft_time[-1]
    time = np.arange(0, 1.0/frequency*N, 1.0/frequency)
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(time, ampl_real)
    plt.grid(True)
    plt.title('Inverse Fast Fourier Transform')
    #plt.xlim([xmin, xmax])
    #plt.ylim([-1, 2.5])
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Inverse_Fast_Fourier_Transform.jpg'
    plt.savefig(file_name)
   
    return time, ampl_real, file_name
 
#==============================================================================
# 実用例
#==============================================================================
# 波形の読込
time, ampl, input_wave_name = ReadWaveCsv()
 
# 高速フーリエ変換
fft_time, fft_ampl, fft_name = fft(time, ampl)
 
# オペレーション・アンプリファイヤー(cut_off 以下の amplitude をゼロにする)
fft_time, fft_ampl, opeamped_fft_name = operation_amplifier(fft_time, fft_ampl, 0.05)
 
# 逆高速フーリエ変換
ret_time, ret_ampl, ifft_name = ifft(fft_time, fft_ampl)
 
# グラフ表示
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 〜 5 に示します。

図1 Input Wave


図2 Fast Fotrier Transform


図3 Operation Amplifier


図4 Inverse Fast Fotrier Transform


図5 Noise Cleared Wave and Original Wave



コメント

この方法の方が微弱(振幅の小さい)な信号に対する処理はやりやすそうです。

ピークの振幅に対して n% 以下はカットするなど使い勝手が良さそうです。

ローパスフィルターも 1 次の周期( 1 番振幅の大きな周期)に対して上下 n% の周期帯から外れた周期をカットするなど、振幅特性と周波数特性についての使い方も何となく分かりました。



以上

Python で信号のノイズ除去 〜ローパスフィルター〜



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


この記事の目次



背景・目的


波形データなどの信号にはノイズがつきものです。

以前の記事で補間や曲線近似などを勉強してきましたが、入手した信号がノイズだらけで補間や近似以前の問題だ…なんて事はよくある事なのではないでしょうか。

そこで、ノイズ除去をやってみたいと思います。

テスト用の簡単な波形データを作り、 その曲線に誤差は与えます。

そして、 ローパスフィルターというフィルターを通して、元の簡単な波形に近づけてみたいと思います。

ローパスフィルターは"ロー"を"パス"するので低い周波数の振幅はそのまま残して高い周波数の振幅をゼロにします。



動作環境


・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 fft(time, ampl):
   
    # サンプリング周期[sec]の計算 #################################################
    sampling_cycle = time[1] - time[0]
   
    # データ数カウント ############################################################
    N = len(ampl)
   
    # 高速フーリエ変換（ FFT ） ####################################################
    fft_ampl = np.fft.fft(ampl)
   
    # FFT の複素数結果を絶対に変換 ###############################################
    abs_fft_amp = np.abs(fft_ampl)
   
    # 振幅をもとの信号に揃える #####################################################
    abs_fft_amp    = abs_fft_amp    / N * 2   # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2       # 直流成分
   
    # 周波数軸のデータ作成 #######################################################
    frequency = np.linspace(0, 1.0/sampling_cycle, N) # 周波数軸　linspace(開始, 終了, 分割数)
   
    # ピーク検出 ################################################################
    maximal_idx = signal.argrelmax(abs_fft_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(frequency[:int(N/2)+1], abs_fft_amp[:int(N/2)+1])
    plt.scatter(frequency[maximal_idx], abs_fft_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 frequency, fft_ampl, file_name
 
#==============================================================================
# ローパスフィルタリング
# cut_off [Hz] 以上の周波数の amplitude をゼロにする
#==============================================================================
def low_pass_filter(fft_time, fft_amp, cut_off):
   
    # データ数カウント ############################################################
    N = len(fft_amp)
   
    # cut_off 以下の周波数の amplitude をゼロにする ################################
    cut_off2 = fft_time[-1] - cut_off
    fft_amp[((fft_time > cut_off)&(fft_time < cut_off2))] = 0 + 0j
          
    # グラフ用に実波形データに変換 #################################################
   
    # FFT の複素数結果を絶対に変換
    abs_fft_amp = np.abs(fft_amp)
   
    # 振幅をもとの信号に揃える
    abs_fft_amp    = abs_fft_amp    / N * 2 # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2     # 直流成分
   
    # ピーク検出 ################################################################
    #maximal_idx = signal.argrelmax(abs_fft_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
    maximal_idx = signal.argrelmax(abs_fft_amp, order=1) # 極大値インデックスの取得
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    #plt.plot(fft_time[:int(N/2)+1], abs_fft_amp[:int(N/2)+1])
    #plt.scatter(fft_time[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.plot(fft_time, abs_fft_amp)
    plt.scatter(fft_time[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.grid(True)
    plt.title('Low Pass Filter')
    plt.xlabel('freqency[Hz]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Low_Pass_Filter.jpg'
    plt.savefig(file_name)
   
    return fft_time, fft_amp, file_name
 
#==============================================================================
# 波形を逆高速フーリエ変換する関数
#==============================================================================
def ifft(fft_time, fft_ampl):
   
    # データ数カウント ############################################################
    N = len(fft_ampl)
   
    # 逆高速フーリエ変換
    ampl = np.fft.ifft(fft_ampl)
   
    # 実数部分のみ
    ampl_real = ampl.real
   
    # 時間軸のデータ作成 ########################################################
    frequency = fft_time[-1]
    time = np.arange(0, 1.0/frequency*N, 1.0/frequency)
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(time, ampl_real)
    plt.grid(True)
    plt.title('Inverse Fast Fourier Transform')
    #plt.xlim([xmin, xmax])
    #plt.ylim([-1, 2.5])
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Inverse_Fast_Fourier_Transform.jpg'
    plt.savefig(file_name)
   
    return time, ampl_real, file_name
 
#==============================================================================
# 実用例
#==============================================================================
# 波形の読込
time, ampl, input_wave_name = ReadWaveCsv()
 
# 高速フーリエ変換
fft_time, fft_ampl, fft_name = fft(time, ampl)
 
# ローパスフィルタリング(任意の周波数以上の amplitude をゼロにする)
fft_time, fft_ampl, passed_fft_name = low_pass_filter(fft_time, fft_ampl, 0.17)
 
# 逆高速フーリエ変換
ret_time, ret_ampl, ifft_name = ifft(fft_time, fft_ampl)
 
# グラフ表示
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 〜 5 に示します。

図1 Input Wave


図2 Fast Fotrier Transform


図3 Low Pass Filter


図4 Inverse Fast Fotrier Transform


図5 Noise Cleared Wave and Original Wave



コメント

意外とキレイに除去されています。

複雑な波形に対して効果があるのか。

まだまだノイズ除去について勉強していきます‼︎



以上

Python で 放射基底関数 (RBF) 補間



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


この記事の目次



背景・目的


前の記事でスプライン曲線で曲線近似を試しました。近似というかほぼ補間ですが…

しかし、放射基底関数 (RBF) の補間を忘れていました。

しかも、 Scipy には放射基底関数 (RBF) の補間が関数として用意されているではありませんか‼︎

この補間方法は絶対記事にして残しておきたい‼︎



動作環境


・Windows 7
・winpython 64bit 3.4.4



プログラム

ソースコード


###############################################################################
# 放射基底関数補間プログラム
###############################################################################
#==============================================================================
# ライブラリインポート
#==============================================================================
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
 
#==============================================================================
# RBF 補間関数
#==============================================================================
def rbf2d(x, y, kernel):
    # 線形 RBF => linear
    # ガウシアン RBF => gaussian
    # 多重二乗 RBF => multiquadric
    # 逆二乗 RBF => inverse
    # 多重調和スプライン RBF => cubic(3次), quintic(5次), thin_plate(薄板スプライン)
   
    # RBF 関数の作成
    rbf = interpolate.Rbf(x, y, function=kernel)
   
    # 入力されている x 値の最小値から最大値の間で任意の間隔で分割する
    xi = np.linspace(x.min(), x.max(), 314)
   
    # xi に対応する yi を計算する
    yi = rbf(xi)
   
    # グラフ描画
    plt.figure(figsize=(10, 8))
    plt.plot( x,  y, 'ro', label='original data point', markersize=4)
    plt.plot(xi, yi,  'r', label='rbf - ' + kernel)
   
    # 格子の設定
    plt.xlim([xi.min(), xi.max()])
    plt.ylim([yi.min()*1.3, yi.max()*1.3])
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
    plt.title("interpolate rbf - " + kernel)
    plt.legend(loc='upper right') # 凡例の位置設定
    # グラフ出力
    file_name = "interpolate rbf - " + kernel + ".jpg"
    plt.savefig(file_name)
    plt.show()
 
#==============================================================================
# 重ね合わせグラフ描画関数
#==============================================================================
def graph(x, y):
    # 線形 RBF => linear
    # ガウシアン RBF => gaussian
    # 多重二乗 RBF => multiquadric
    # 逆二乗 RBF => inverse
    # 多重調和スプライン RBF => cubic(3次), quintic(5次), thin_plate(薄板スプライン)
   
    # カーネルの種類配列作成
    kernel = np.array(['linear', 'gaussian', 'multiquadric', 'inverse', 'cubic', 'quintic', 'thin_plate'])
   
    # グラフの作成
    plt.figure(figsize=(10, 8))
   
    for i in range(len(kernel)):
        # RBF 関数の作成
        rbf = interpolate.Rbf(x, y, function=kernel[i])
       
        # 入力されている x 値の最小値から最大値の間で任意の間隔で分割する
        xi = np.linspace(x.min(), x.max(), 314)
       
        # xi に対応する yi を計算する
        yi = rbf(xi)
       
        # グラフ描画
        plt.plot(xi, yi, label='rbf - ' + kernel[i])
   
    # 格子の設定
    plt.grid(which='major', color='black', linestyle='-')
    plt.grid(which='minor', color='black', linestyle='-')
    # 格子の設定
    plt.xlim([xi.min(), xi.max()])
    plt.ylim([yi.min()*1.4, yi.max()*1.4])
    plt.title("interpolate rbf")
    plt.legend(loc='upper right') # 凡例の位置設定
    # グラフ出力
    file_name = "interpolate rbf.jpg"
    plt.savefig(file_name)
    plt.show()
 
#==============================================================================
# テスト
#==============================================================================
# 値の準備 =====================================================================
N = 2**5                           # サンプリング数
dt = 0.01                          # サンプリング周期
ampl = np.array([  2])             # 振幅
freq = np.array([  1])             # 周波数
peri_axis = np.arange(0, N*dt, dt) # 周期軸
error = 0.01                       # 誤差係数
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) # 誤差を与えた入力波形
 
# 補間 ＆ 出力 =================================================================
rbf2d(peri_axis, y, 'linear')
rbf2d(peri_axis, y, 'gaussian')
rbf2d(peri_axis, y, 'multiquadric')
rbf2d(peri_axis, y, 'inverse')
rbf2d(peri_axis, y, 'cubic')
rbf2d(peri_axis, y, 'quintic')
rbf2d(peri_axis, y, 'thin_plate')
graph(peri_axis, y)



結果

放射基底関数 (RBF) 補間の結果を図1 〜 図7 に示します。

また、全てを重ね合わせた結果を図8 に示します。

図1 線形 RBF


図2 ガウシアン RBF


図3 多重二乗 RBF


図4 逆二乗 RBF


図5 多重調和スプライン RBF(3次)


図6 多重調和スプライン RBF(5次)


図7 多重調和スプライン RBF(薄板)


図8 重ね合わせ



コメント

補間結果ですが、若干振動してますね。
ガウシアンがわりと大きく振動しててなんかショックでした。

自動でノイズ除去をして、自動でスプラインや放射基底関数 (RBF) 補間を選択して、自動で補間してキレイな波形データが得られる様になったら、良いんですけどね…

曲面の補間もやってみたいですね。



以上

Python で曲線近似（フィッティング）テスト用データの作成



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


この記事の目次



背景・目的


曲線近似のプログラムを作っていく上で、テストは欠かせません。

もちろん、どんなプログラムを作るにしてもテストは必要ですが…

そこで、テスト用のある程度複雑な曲線を作り、 その曲線上の点をランダムに抽出して、誤差を与えます。

そして、 Microsoft Word ファイルと csv ファイルにそれぞれ曲線仕様、抽出点仕様および抽出点の座標データを出力するプログラムを作ります。



動作環境


・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_graph(x, y, x1, y1, xmin, xmax, ymin, ymax):
    plt.figure(figsize=(5.5, 3.5))
   
    # 元々の値の出力
    plt.plot(x, y, 'ro', label="original data point")
    plt.plot(x, y, 'r')
   
    # 誤差を含んだ値の出力
    plt.plot(x1, y1, 'bo', 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=(5.5, 3.5))
   
    # 各 x 値における誤差の出力
    plt.bar(x, y, linewidth = 0, width = 0.0007)
    plt.plot(x, y, 'ro', 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
 
#==============================================================================
# テスト用データを作成する関数 Wave の実行結果報告書作成関数
#==============================================================================
def DocumentReport(wave_number, data_number, dt, ampl, freq, t, f, img_name_wave, img_name_err):
   
    # add_heading の第２引数
    # 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に示す。また、付加した誤差を図2に示す。なお、テスト用データのデジタルデータをまとめる。')
   
    doc.add_picture(img_name_wave) # 作成したテスト用データグラフの挿入
    p = doc.add_paragraph('図1　テスト用データ\n')
    p.alignment = WD_ALIGN_PARAGRAPH.CENTER
   
    doc.add_picture(img_name_err) # 付加した誤差グラフの挿入
    p = doc.add_paragraph('図2　付加した誤差\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\\wave for test.docx')
 
#==============================================================================
# 作成したテスト用データのデジタルデータを csv 出力する関数
#==============================================================================
def WriteWaveCsv(t, f):
    w = np.stack([t, f])
    np.savetxt('C:\\Users\\UserName\\Desktop\\wave for test.csv', w.T, delimiter=',')
 
#==============================================================================
# テスト用データ作成
#==============================================================================
# 値の準備 =====================================================================
N = 2**5                           # サンプリング数
dt = 0.001                         # サンプリング周期
ampl = np.array([ 1,  1,  2])      # 振幅
freq = np.array([60, 45, 100])     # 周波数
peri_axis = np.arange(0, N*dt, dt) # 周期軸
error = 0.1                        # 誤差係数
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) # グラフの ｘ 軸の最小値設定
xmax = max(peri_axis) + (3*dt) # グラフの ｘ 軸の最大値設定
ymin = min(y) * 1.5            # グラフの y 軸の最小値設定
ymax = max(y) * 1.5            # グラフの 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
 
# グラフ出力 ====================================================================
img_name_wave = wave_graph(peri_axis, origin_y, peri_axis, y, xmin, xmax, ymin, ymax)
img_name_err  = err_graph(peri_axis, err, xmin, xmax, ymin_err, ymax_err)
 
# テスト用データ作成報告書の作成 ===================================================
DocumentReport(len(ampl), N, dt, ampl, freq, peri_axis, y, img_name_wave, img_name_err)
 
# テスト用データの csv 出力 =======================================================
WriteWaveCsv(peri_axis, y)



結果

このプログラムから出力される曲線近似（フィッティング）テスト用データのグラフ、テスト用データに付加した誤差のグラフ、 Word ファイル、 csv ファイルをそれぞれ図1 〜 図4 に示します。

図1 テスト用データグラフ


図2 誤差グラフ


図3 Word ファイル（冒頭）


図4 csv ファイル（冒頭）



コメント

曲線近似プログラムをテストするときはこのプログラムでテスト用データを大量に作成して、即座にテストに着手できるようになったかと思います。

以前の記事で波形のデータを作成するプログラムを作ったことがありました。

しかし、作成した波形データに誤差を与え忘れている事に気がつきました。



以上

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 の第２引数
    # 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 の精度向上について調べていきたいです。



以上