Tremor analysis notebook.ipynb
This notebook helps analyze tremor data via a fast fourier transform, to report the frequency and magnitude of a tremor.
This notebook expects to analyze a data file from the "Physics Toolbox Sensor Suite" phone app.
target_file
: CSV data file name.Before doing a fourier transform, data should be a uniform sampling – there should be an equal amount of time between each data timepoint. This isn't always true in the raw data. Use an interpolation to get uniform data.
This notebook expects to analyze a data file from the "Physics Toolbox Sensor Suite" phone app.
target_file
: CSV data file name.target_file = "sensor-4.csv"
# Edit this to change the plot title. The default uses your file name.
plot_title = "FFT tremor analysis for {}".format(target_file)
import math
import statistics
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.interpolate import interp1d
import seaborn
data_raw = pd.read_csv(target_file, parse_dates=['time'])
Before doing a fourier transform, data should be a uniform sampling – there should be an equal amount of time between each data timepoint. This isn't always true in the raw data. Use an interpolation to get uniform data.
# Get the median interval time for this data (i.e. typical interval)
intervals = [data_raw.index[i] - data_raw.index[i-1] for
i in range(1, len(data_raw.index))]
interval = statistics.median(intervals)
# Convert to unix epoch timestamp for easier math.
ts_raw = [t.timestamp() for t in data_raw.time]
# Set up interpolation for total gForce.
gf_raw = list(data_raw.gFTotal)
interpolate = interp1d(ts_raw, gf_raw)
# Create uniform timepoints, derive interpolated gForce values.
ts_uniform = np.linspace(ts_raw[0], ts_raw[-1], len(ts_raw))
avg_delta = (ts_raw[-1] - ts_raw[0]) / len(ts_raw)
gf_uniform = interpolate(ts_uniform)
gf_fft_all = np.fft.fft(gf_uniform)
freqs_all = np.fft.fftfreq(len(ts_uniform), d=avg_delta)
# discard complex conjugate
target_len = int(len(freqs_all)/2)
freqs = freqs_all[1:target_len]
gf_fft = gf_fft_all[1:target_len]
plt.figure()
plt.plot(freqs, gf_fft)
plt.title(plot_title)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
# Get the maximum value, report the frequency and magintude.
peak_index = np.argmax(gf_fft)
peak_freq = freqs[peak_index]
peak_magnitude = abs(gf_fft[peak_index])
print("Peak frequency: {} Hz".format(peak_freq))
print("Peak magnitude: {}".format(peak_magnitude))