EEG Hacker, 2015-01-17, MIT License
OpenBCI data file. Do some simple visualizations to make sure that it looks OK.
First, we have to import some libraries...though this might be done automagically for IPython Notebooks...I'm not sure.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
# Define Data to Load...be sure to UNZIP the data!!!
pname = 'SavedData/'
test_id = 0;
if (1):
#day 2...after caffeine
fname = 'OpenBCI-RAW-2015-01-25_18-28-58_eyesOpen_AMSine.txt';
t_lim_sec = [110.0, 280.0]
test_id = 1
elif (0):
#day 2...after caffeine
fname = 'OpenBCI-RAW-2015-01-25_19-04-27_twoTone_focus2500_focus1000.txt'; #2500 is 43 Hz, 1000 is 37 Hz
t_lim_sec = [42.0, 92.0]
test_id = 2
else:
#day 2...after caffeine
fname = 'OpenBCI-RAW-2015-01-25_19-15-48_twoTone_eyesClosed_focus1000_focus2500.txt';
t_lim_sec = [85.0, 135.0]
test_id = 3
# 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))
Number of discontinuities in the packet counter: 13
# 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(eeg_data_uV[:, 0])) / fs_Hz
# plot each channel of EEG data...raw
nchan = 8 #normally 8 or 16
ncol = 2;
nrow = nchan / ncol
plt.figure(figsize=(ncol*8, nrow*4.5))
for Ichan in range(nchan):
plt.subplot(nrow,ncol,Ichan+1)
plt.plot(t_sec,eeg_data_uV[:, Ichan])
plt.xlabel("Time (sec)")
plt.ylabel("Raw EEG (uV)")
plt.title("Channel " + str(Ichan+1))
#plt.ylim([-1.5, 1.5])
plt.tight_layout()
NFFT = 256*2 # pitck the length of the fft
FFTstep = 0.5*fs_Hz # do a new FFT every half second
overlap = NFFT - FFTstep # half-second steps
f_lim_Hz = [0, 50] # frequency limits for plotting
plt.figure(figsize=(ncol*8, nrow*4.5))
for Ichan in range(nchan):
ax = plt.subplot(nrow,ncol,Ichan+1)
data = np.array(eeg_data_uV[:,Ichan])
data = data - np.mean(data,0)
spec_PSDperHz, freqs, t = mlab.specgram(data,
NFFT=NFFT,
window=mlab.window_hanning,
Fs=fs_Hz,
noverlap=overlap
) # returns PSD power per Hz
spec_PSDperBin = spec_PSDperHz * fs_Hz / float(NFFT) #convert to "per bin"
plt.pcolor(t, freqs, 10*np.log10(spec_PSDperBin)) # dB re: 1 uV
#plt.clim(20+np.array([-45, 0]))
plt.clim(20-7.5-3.0+np.array([-30, 0]))
#plt.xlim(t_sec[0], t_sec[-1])
plt.xlim(np.array(t_lim_sec)+np.array([-10, 10]))
#plt.ylim([0, fs_Hz/2.0]) # show the full frequency content of the signal
plt.ylim(f_lim_Hz)
plt.xlabel('Time (sec)')
plt.ylabel('Frequency (Hz)')
plt.title("Channel " + str(Ichan+1))
# add annotation for FFT Parameters
ax.text(0.025, 0.95,
"NFFT = " + str(NFFT) + "\nfs = " + str(int(fs_Hz)) + " Hz",
transform=ax.transAxes,
verticalalignment='top',
horizontalalignment='left',
backgroundcolor='w')
plt.colorbar()
plt.tight_layout()
#do the channel of interest
Ichan = 7-1
#get the time period of interest
if (test_id == 1):
if (0):
t_focus_sec = 182.0+60.0+10+np.array([0, 20])
elif (1):
#should be 40 Hz?
t_focus_sec = np.array([205.0, 220.0]) #see in chans 3, 6, and 7?
elif (0):
#at both 37 and 43 Hz
t_focus_sec = 250.0+np.array([0, 20]) #see in chan 6?
elif (test_id == 2):
if (1):
#both 37 and 43, stronger on 43?
t_focus_sec = 43.0+np.array([0, 18])
else:
#both 37 and 43, stronger on 37?
t_focus_sec = 43.0+30.0+np.array([0, 18])
elif (test_id == 3):
if (1):
#both 37 and 43, stronger 37?
t_focus_sec = 87.0+np.array([0, 15])
else:
#both 37 and 43, stronger 43?
t_focus_sec = 115.0+np.array([0, 15])
for Ichan in range(nchan):
#get the time period of interest
t_focus_samp = t_focus_sec * fs_Hz + 1
data = np.array(eeg_data_uV[t_focus_samp[0]:t_focus_samp[1],Ichan])
#convert to the frequency domain
NFFT = 256*2
data = data - np.mean(data,0)
spec_PSDperHz, freqs, t = mlab.specgram(data,
NFFT=NFFT,
window=mlab.window_hanning,
Fs=fs_Hz,
noverlap=overlap
) # returns PSD power per Hz
spec_PSDperBin = spec_PSDperHz * fs_Hz / float(NFFT) #convert to "per bin"
#get the mean spectrum in that time
spectrum_PSDperHz = np.mean(spec_PSDperHz,1) #time is horizontal in the 2D array?
#plot
if (Ichan == 0):
plt.figure(figsize=(10, 5))
#plt.figure(figsize=(ncol*10, nrow*5))
ax = plt.subplot(1,1,1)
plt.plot(freqs, 10*np.log10(spectrum_PSDperHz)) # dB re: 1 uV
#plt.xlim([0, fs_Hz/2.0]) # show the full frequency content of the signal
plt.xlim([30, 50])
plt.ylim([-20.0, 5.0])
plt.plot(37.0*np.array([1, 1]),ax.get_ylim(),'k--',linewidth=2)
plt.plot(40.0*np.array([1, 1]),ax.get_ylim(),'k--',linewidth=2)
plt.plot(43.0*np.array([1, 1]),ax.get_ylim(),'k--',linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD per Hz (dB re: 1uV^2/Hz)')
#plt.title("Channel " + str(Ichan+1))
plt.title(fname + ': t = [' + str(t_focus_sec[0]) + ' ' + str(t_focus_sec[1]) + '] sec')
plt.legend(['1', '2', '3', '4', '5', '6', '7', '8'])
# add annotation for FFT Parameters
ax.text(0.025, 0.05,
"NFFT = " + str(NFFT) + "\nfs = " + str(int(fs_Hz)) + " Hz",
transform=ax.transAxes,
verticalalignment='bottom',
horizontalalignment='left',
backgroundcolor='w')
<matplotlib.text.Text at 0xc426c18>
eeg_data_uV.shape
(72443L, 8L)