import numpy as np
from scipy import fftpack, signal
%matplotlib notebook
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
data = np.genfromtxt('data/TGV_data.csv.bz2', delimiter=',', names=True)
time = data['Time_s']
omega_y = data['Gyroscope_y_rads']
_ = plt.plot(time, omega_y)
Before analyzing the data any further, let us take a look at the time elapsed between subsequent measurements.
time_intervals = time[1:]-time[:-1]
_ = plt.hist(time_intervals, bins=100)
We will ignore these small differences and assume the data to be equally spaced in time with the mean time difference between subsequent data points.
delta_t = np.mean(time_intervals)
delta_t
Fourier transformation of a discrete real signal with rfft
:
$$y_j = \sum_{k=0}^{n-1} x_k \exp\left(-\text{i}\frac{2\pi jk}{n}\right)$$
With exception of $y_0$ and, for even $n$, $y_{n/2}$ all $y_j$ are complex.
rfft
returns an array $y_0, \text{Re}(y_1), \text{Im}(y_1),\ldots$
fft_y = fftpack.rfft(omega_y)
_ = plt.plot(fft_y[::2])
_ = plt.plot(fft_y[1::2])
freqs = fftpack.rfftfreq(omega_y.shape[0], delta_t)
_ = plt.plot(freqs[::2], fft_y[::2])
_ = plt.plot(freqs[1::2], fft_y[1::2])
fft_y = np.insert(fft_y, 0, 0)
fft_square = fft_y*fft_y
power = fft_square[::2]+fft_square[1::2]
power = 2*delta_t/fft_square.shape[0]*power
_ = plt.plot(freqs[::2], power)
_ = plt.plot(*signal.periodogram(omega_y, 1/delta_t))
We are interested in the signal around 1.4 Hz. Filter out frequencies beyond 3 Hz.
filter_coeffs = signal.firwin(301, 3, pass_zero=True, fs=1/delta_t)
_ = plt.plot(filter_coeffs, '.')
freqs, response = signal.freqz(filter_coeffs, fs=1/delta_t)
_ = plt.plot(freqs, 20*np.log10(np.abs(response)))
_ = plt.polar(np.angle(response), np.abs(response))
omega_y_filtered = signal.convolve(omega_y, filter_coeffs, mode='valid')
omega_y.shape, omega_y_filtered.shape
_ = plt.plot(time[150:-150], omega_y[150:-150])
_ = plt.plot(time[150:-150], omega_y_filtered)
Filter out anything but the range from 53-70 Hz
filter_coeffs = signal.firwin(301, [53, 70], pass_zero='bandpass', fs=1/delta_t)
freqs, response = signal.freqz(filter_coeffs, fs=1/delta_t)
_ = plt.plot(freqs, 20*np.log10(abs(response)))
omega_y_filtered = signal.convolve(omega_y, filter_coeffs, mode='valid')
_ = plt.plot(time[150:-150], omega_y[150:-150])
_ = plt.plot(time[150:-150], omega_y_filtered)
_ = plt.plot(fftpack.rfftfreq(omega_y_filtered.shape[0], delta_t)[::2],
fftpack.rfft(omega_y_filtered)[::2])
freq, sp_time, Sxx = signal.spectrogram(omega_y, fs=1/delta_t)
_ = plt.pcolormesh(sp_time, freq, Sxx)
fig = plt.figure()
ax = Axes3D(fig)
_ = ax.plot_surface(sp_time, freq[:, np.newaxis], Sxx, cmap=cm.viridis)