integrate qu8k accelerometer data
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import scipy as sp | |
from scipy import fft, ifft | |
from scipy.integrate import cumtrapz | |
from scipy.integrate import simps | |
import csv | |
from matplotlib import rc | |
fig_width_pt = 469.75 | |
inch_per_pt = 1.0 / 72.27 | |
figw = fig_width_pt * inch_per_pt | |
figh = (620.43 / 2.1) * inch_per_pt | |
fig_size = [figw, figh] | |
rc('figure', figsize=fig_size) | |
import matplotlib.pylab as plt | |
def fft_diff(y, L): | |
i = 1.0j | |
n = len(y) | |
k = sp.append(sp.array(range(0, n/2)), 0.0) | |
k = (2.0 * sp.pi / L) * sp.append(k , sp.array(range(-n/2+1,0))) | |
return(ifft(i * k * fft(y)).real) | |
def fft_int(y, x, L): | |
i = 1.0j | |
n = len(y) | |
k = sp.append(sp.array(range(0, n/2)), 0.0) | |
k = (2.0 * sp.pi / L) * sp.append(k , sp.array(range(-n/2+1,0))) | |
fy = fft(y) | |
y_int = -i * fy / k | |
y_int[0] = 0.0 | |
y_int[n/2] = 0.0 | |
return(ifft(y_int).real + (x*fy[0]/n).real) | |
def c_T(T): | |
gamma = 1.4 | |
R = 1716.59 | |
return(sp.sqrt(gamma*R*T)) | |
def c_alt(h): | |
T0 = 518.67 # R | |
theta = sp.ones(h.shape) | |
theta[h<278385.8] = 1.237723 - h[h<278385.8] / 472687.0 | |
theta[h<232939.6] = 1.434843 - h[h<232939.6] / 337634.0 | |
theta[h<167322.8] = 0.939268 | |
theta[h<154199.5] = 0.482561 + h[h<154199.5] / 337634.0 | |
theta[h<104986.9] = 0.682457 + h[h<104986.9] / 945374.0 | |
theta[h<65616.8] = 0.751865 | |
theta[h<36089.2] = 1.0 - h[h<36089.2] / 145442.0 | |
return(c_T(theta*T0)) | |
reader = csv.reader(open("qu8k_ascent.csv","r")) | |
skiplines = 2 | |
header = [] | |
data = [] | |
for i,row in enumerate(reader): | |
if(i<skiplines): | |
header.append(row) | |
else: | |
data.append([float(col) for col in row[0:len(row)-1]]) | |
traj = sp.array(data) | |
traj = traj[0:traj.shape[0]-7,:] | |
traj[:,0] = traj[:,0] + 0.04 | |
# acceleration in the frequency domain: | |
a_f = fft(traj[:,2]) | |
v_fft = fft_int(traj[:,2], traj[:,0], max(traj[:,0])-min(traj[:,0])) | |
# set the initial condition to v=0, arbitrary constant of integration | |
v_fft = v_fft - v_fft[0] | |
d_fft = fft_int(v_fft, traj[:,0], max(traj[:,0])-min(traj[:,0])) | |
d_fft = d_fft - d_fft[0] + 3900.0 # black rock is at 3900 ft | |
mach_fft = v_fft / c_alt(d_fft) | |
v_trapz = sp.append(0.0, cumtrapz(traj[:,2], traj[:,0])) | |
fig = plt.figure() | |
plt.title(r'"Qu8k" at Black Rock, NV 30 Sept 2011') | |
ax1 = fig.add_subplot(111) | |
ax1.plot(traj[:,0], traj[:,2]/1e3, 'r-', label="acceleration") | |
#ax1.plot(traj[:,0], traj[:,3], label="velocity, rectangular") | |
#ax1.plot(traj[:,0], v_trapz, label="velocity, trapezoidal") | |
ax1.plot(traj[:,0], v_fft/1e3, 'g-', label="velocity") | |
ax1.plot(traj[:,0], mach_fft, 'm-', label="Mach") | |
xmin,xmax,ymin,ymax = ax1.axis() | |
ax1.set_xlim([0,xmax]) | |
ax1.legend(loc=9) | |
ax1.set_xlabel(r'time $(s)$') | |
ax1.set_ylabel(r'$kft/s$ or $kft/s/s$ or $M$') | |
ax2 = ax1.twinx() | |
ax2.plot(traj[:,0], d_fft/1e3, 'b-', label='altitude') | |
ax2.set_ylabel(r'altitude $(kft)$', color='b') | |
ax1.text(20,0.1,'accelerometer data:\n ddeville.com/derek/Qu8k.html') | |
plt.savefig("qu8k_trajectory.pdf") | |
plt.figure() | |
plt.title(r'difference from FFT-based integration') | |
plt.semilogy(traj[:,0], abs(traj[:,3]-v_fft), label=r'1st order (rectangular)') | |
plt.semilogy(traj[:,0], abs(v_trapz-v_fft), label=r'2nd order (trapezoidal)') | |
plt.legend(loc=0) | |
plt.ylabel(r'velocity magnitude $(ft/s)$') | |
plt.xlabel(r'time $(s)$') | |
plt.savefig("qu8k_trajectory_error.pdf") | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment