# First, import the Python libraries that we need
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# Define Data to Load...be sure to UNZIP the data!!!
pname = 'SavedData/'
fname = 'OpenBCI-RAW-2014-11-23_18-54-57.txt'

# load data into numpy array
fs_Hz = 250.0        # assumed sample rate for the EEG data
data = np.loadtxt(pname + fname,
                  delimiter=',',
                  skiprows=5)

# check the packet counter for dropped packets
data_indices = data[:, 0]   # the first column is the packet index
d_indices = data_indices[2:]-data_indices[1:-1]
n_jump = np.count_nonzero((d_indices != 1) & (d_indices != -255))
print("Number of discontinuities in the packet counter: " + str(n_jump))

# parse the data out of the values read from the text file
# eeg_data_uV = data[:, 1:(8+1)] # EEG data, microvolts
accel_data_counts = data[:, 9:(11+1)]  # note, accel data is NOT in engineering units

# convert units
unused_bits = 4 #the 4 LSB are unused?
accel_data_counts = accel_data_counts / (2**unused_bits) # strip off this extra LSB
scale_G_per_count = 0.002  # for full-scale = +/- 4G, scale is 2 mG per count
accel_data_G = scale_G_per_count*accel_data_counts # convert to G

# create a vector with the time of each sample
t_sec = np.arange(len(accel_data_G[:, 0])) / fs_Hz

# The accel data seems to be padded with lots of zeros. 
I = np.nonzero(accel_data_G[:,0])
I=I[0]
print("Accel data exists in " + str(len(I)) + 
      " of " + str(len(t_sec)) + \
      " (" + str((100*len(I))/len(t_sec)) + 
      "%) of data packets");

# Only get the non-zero entries
accel_data_G = accel_data_G[I,:]
t_sec = t_sec[I]

# plot each channel of accelerometer data
for Iaxis in range(3):      
    plt.figure(figsize=(9.0, 2.5))
    plt.plot(t_sec,accel_data_G[:, Iaxis])
    plt.plot([t_sec[0], t_sec[-1]],[0, 0],'k--')
    plt.plot([t_sec[0], t_sec[-1]],[1.0, 1.0],'k:')
    plt.plot([t_sec[0], t_sec[-1]],[-1.0, -1.0],'k:')
    plt.xlabel("Time (sec)")       
    plt.ylabel("Acceleration (G)")
    plt.title("Accelerometer Channel " + str(Iaxis))
    plt.ylim([-1.5, 1.5])
    
    # add annotations to show each movement
    if (Iaxis==0) :
        plt.text(1.0,0.25,"Flat and Level",
                 horizontalalignment='left',
                 backgroundcolor='w',
                 verticalalignment='center')
        plt.text(10.5,-0.8,"Roll Right",
                 horizontalalignment='right',
                 backgroundcolor='w',
                 verticalalignment='center')
        plt.text(18.5,0.8,"Roll Left",
                 horizontalalignment='left',
                 backgroundcolor='w',
                 verticalalignment='center')
    elif (Iaxis==1) :
        plt.text(20,-0.8,"Nose Down",
                 horizontalalignment='right',
                 backgroundcolor='w',
                 verticalalignment='center')
        plt.text(27,0.8,"Nose Up",
                 horizontalalignment='left',
                 backgroundcolor='w',
                 verticalalignment='center')
    elif (Iaxis==2) :
        plt.text(30.5,-0.8,"Upside Down",
                 horizontalalignment='right',
                 backgroundcolor='w',
                 verticalalignment='center')
        

# compute magnitude of the acceleration vector
mag_accel_G = np.sqrt(accel_data_G[:,0]**2 + 
                      accel_data_G[:,1]**2 + 
                      accel_data_G[:,2]**2)
median_mag_accel_G = np.median(mag_accel_G)

#plot magnitude of acceleration
plt.figure(figsize=(9, 2.5))
plt.plot(t_sec,mag_accel_G)
plt.plot([t_sec[0], t_sec[-1]],[1.0, 1.0],'k--')
plt.xlabel("Time (sec)")       
plt.ylabel("Acceleration (G)")
plt.title("Magnitude of Acceleration")
plt.ylim([0.6, 1.4])
plt.text(1,1.3,"Median = " + "{:.3f}".format(median_mag_accel_G) + "G")
plt.show()

print "Median acceleration should be near 1.0 G"
print "Median acceleration is actually " + \
       "{:.3f}".format(median_mag_accel_G) + " G"