Python

【Raspberry Pi】加速度センサーの時間信号をフーリエ変換(FFT)で周波数表示

rasbperry PiでFFT

こんにちは！　けい(@Keimameshiba)です．

今回は，加速度センサーから時間信号を取得して，FFTにより周波数成分で表示する方法についてまとめていきます．

前回の 【RASPBERRY PI】加速度センサーの値をCSVファイルに出力して可視化(MPU6050) の続きなので，まだ読んでない方はぜひ読んでみてください．

【Raspberry Pi】加速度センサーの値をCSVファイルに出力して可視化(MPU6050)

ソースコード

# -*- coding: utf-8 -*-

import smbus
import math
from time import sleep
import time
import numpy as np

DEV_ADDR = 0x68

ACCEL_XOUT = 0x3b
ACCEL_YOUT = 0x3d
ACCEL_ZOUT = 0x3f
TEMP_OUT = 0x41
GYRO_XOUT = 0x43
GYRO_YOUT = 0x45
GYRO_ZOUT = 0x47

PWR_MGMT_1 = 0x6b
PWR_MGMT_2 = 0x6c   

bus = smbus.SMBus(1)
time.sleep(0.2)
# 初期設定
bus.write_byte_data(DEV_ADDR, 0x6B, 0x80) # RESET
time.sleep(0.25)
bus.write_byte_data(DEV_ADDR, 0x6B, 0x00) # RESET
time.sleep(0.25)
bus.write_byte_data(DEV_ADDR, 0x6A, 0x07) # RESET
time.sleep(0.25)
bus.write_byte_data(DEV_ADDR, 0x6A, 0x00) # RESET
time.sleep(0.25)
bus.write_byte_data(DEV_ADDR, 0x1A, 0x00) # CONFIG
bus.write_byte_data(DEV_ADDR, 0x1B, 0x18) # +-2000/s
bus.write_byte_data(DEV_ADDR, 0x1C, 0x08) # +-4g


def read_word(adr):
    high = bus.read_byte_data(DEV_ADDR, adr)
    low = bus.read_byte_data(DEV_ADDR, adr+1)
    val = (high << 8) + low
    return val

def read_word_sensor(adr):
    val = read_word(adr)
    if (val >= 0x8000):  return -((65535 - val) + 1)
    else:  return val

def get_temp():
    temp = read_word_sensor(TEMP_OUT)
    x = temp / 340 + 36.53      # data sheet(register map)記載の計算式.
    return x

def getGyro():
    x = read_word_sensor(GYRO_XOUT)/ 16.4
    y = read_word_sensor(GYRO_YOUT)/ 16.4
    z = read_word_sensor(GYRO_ZOUT)/ 16.4
    return [x, y, z]


def getAccel():
    x = read_word_sensor(ACCEL_XOUT)/ 8192
    y= read_word_sensor(ACCEL_YOUT)/ 8192
    z= read_word_sensor(ACCEL_ZOUT)/ 8192
    return [x, y, z]

# data parameta
N = 512
dt = 0.02                     
t = np.arange(0, N*dt, dt)
freq = np.linspace(0, 1.0/dt, N)

time_data = np.array([])
count = 0
while(True):
    for i in range(N):
        ax, ay, az = getAccel()
        gx, gy, gz = getGyro()
        print(count)
        count += 1
        time_data = np.append(time_data, az)
        time.sleep(0.02)
    
    #FFT
    F = np.fft.fft(time_data)
    Amp = np.abs(F)
    break

file = open("accel_data_fft_time.csv", "w")
for i,val in enumerate(time_data):
    file.write(str(i*dt) + "," + str(val) + "\n")
file.close()

fft_file = open("accel_data_fft_fre.csv", "w")
for i,val in enumerate(Amp):
    fft_file.write(str(i/(dt*N)) + "," + str(val) + "\n")
fft_file.close()

時間信号の取得

ここは前回と全く同じなので，気になる方は前回の記事を読んでください．

ax, ay, az = getAccel()
gx, gy, gz = getGyro()

ここの２行で，加速度とジャイロをセンサーから読み取っています．

今回のコードでは，azつまりZ方向の加速度のみを用います．

73行目のtime_dataに，加速度の時間信号を追加していきます．

高速フーリエ変換【FFT】

信号処理の鉄板手法のFFTを行います．FFTは離散フーリエ変換(DFT)の計算を，高速に行うための手法のことで，DFTと原理は同じです．(気になる人はこの記事がおすすめ)

FFTを簡単に行えるライブラリが，numpyにあるのでそちらを用います．FFTには各種パラメータが必要なので，それを定義しているのが67行目からです．

順に，

  1. 時間信号が何個のデータあるか：N
  2. サンプリング周期：dt
  3. 別途，グラフ表示したい時の時間配列：t
  4. 別途，グラフ表示したい時の周波数配列：freq

今回はtとfreqは使っていませんが，matplotlibなどで時間信号と周波数成分を表示したい時に，それぞれが第１軸として表示できます．

np.fft.fft()関数に，一次元配列を渡すことでフーリエ変換できます．返り値は複素数なので，絶対値を取ってAmpに格納します．

CSVファイルに保存して，エクセルで可視化

加速度の値の，時間成分と周波数成分が作れたので，それぞれをcsvファイルに保存します．

時間信号のファイルをaccel_data_fft_time.csvというファイル名，周波数成分のファイルをaccel_data_fft_fre.csvというファイル名で保存します．

forの繰り返しに，enumerate関数を組み合わせてループを回します．enumerate は一つ目のiにインデックス，二つ目のvalにenumerate関数の中身の値を返します．

なので例えば，time_dataに[25, 30, 35, 40]という配列が格納されていたら，forループでiには順に0,1,2,3，valには順に25,30,35,40が格納されます．

また，エクセルでそのままｘ軸とｙ軸で表示したいので，インデックスに１つのデータ幅をかけてｘ軸の値にしています．

こちらがcsvにファイルに保存して，そのままグラフで表示したものです．

１秒間に１回の間隔で振動させたのですが，保存したデータの10秒間に13個の振動が測定されました．これは実際には20ms間隔で測定できず，より多くの時間がかかってしまったため，より多くの振動が測定された，ということだと思います．

なぜなら，time.sleep関数を用いて20ms待つことで，サンプリング周期を20msにしたからです．これでは，内部の計算時間を0sと近似しているので遅れが生じました．

まとめ

今回は，加速度センサーの時間信号を周波数成分で表示してみました．

より正確に行うために，サンプリング周期をより正確にする必要があることが分かりました．

FFTを活用すれば，色々なことに応用できるので知っていて損はないと思います．