以前、ESP8266と3軸加速度センサーADXL345を使って振動を測定し、Ambientで可視化しました。

今回は、このデーターを使い、振動データーの周波数成分を調べます。

この振動データーは3軸加速度センサーADXL345をテーブルの上に置き、横2方向(x軸、y軸)と垂直方向(z軸)に机を軽く叩き、x軸、y軸、z軸、3方向の振動を10m秒間隔(100Hz)で10秒間サンプリングしたものです。データーは以下のチャネルで公開しています。

このデーターを高速フーリエ変換(FFT)処理して、周波数成分を調べます。今回はPythonのFFTライブラリーを使うことにしました。まず、必要なライブラリーをインポートします。

import numpy as np
import scipy.fftpack
from pylab import *
import pandas as pd

Ambientに蓄積したデーターをPythonで読み込むには、Ambientライブラリーを使う方法と、Ambientのデーターをcsvファイルにダウンロードし、Pythonでcsvファイルを読み込む方法があります。今回はcsvファイルから読み込んでみます。

df = pd.read_csv('data/data1475.csv')

こんな感じでpandasのDataFrameとして読み込めます。中身を確認します。

df.head(3)

ここから、z軸加速度のデーターを抜き出し、最大値、最小値を求め、値が-1から1の範囲になるようの正規化します。

x = df.iloc[0:, 3:4].values.flatten()
(np.max(x), np.min(x)) # 結果は(390.0, -253.5)
x = [y / 390.0 for y in x]

データーを確認します。

plot(x[0:1000])
show()

この中から616サンプル目から60サンプル、時間にすると測定開始6.16秒目から0.6秒の振動データーの周波数成分をFFT処理してみます。

 

start = 616 # サンプリングする開始位置 
N = 60 # FFTのサンプル数 
fs = 100.0 # サンプリング周波数 100Hz => 10m秒間隔でサンプリング 

hanningWindow = np.hanning(N) # ハニング窓 
X = scipy.fftpack.fft(hanningWindow * x[start:start+N]) 

freqList = scipy.fftpack.fftfreq(N, d=1.0/ fs) # 周波数軸の値を計算 

amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X] # 振幅スペクトル

元の波形と振幅スペクトルを表示します。


# 元の波形 
plot(range(start, start+N), x[start:start+N]) 
axis([start, start+N, -1.0, 1.0]) 
xlabel("time [sample]") 
ylabel("amplitude") 
show() 

# 振幅スペクトル 
plot(freqList, amplitudeSpectrum, marker= 'o', linestyle='-') 
axis([0, fs/2, 0, 1]) 
xlabel("frequency [Hz]") 
ylabel("amplitude spectrum") 
show()


30Hzあたりにピークのある周波数成分の振動であることが確認できました。

元のデーターは10m秒周期で10秒間サンプリングしています。これを先頭から128サンプルFFT処理し、16サンプルずつずらしながら、最後までFFT処理し、周波数成分の変化を見ると、次のようになりました。

振動データーを継続的に測定し、周波数成分の変化を抽出することで、工場設備の故障予知などにつなげていけそうです。

FFT処理は次のサイトを参考にさせていただきました。