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

この記事の目次

背景・目的

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" 関数の中でかなりテキトーにやっているので、信頼できない類似度です。

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