created flag in feature engineer to concatenate imaginary values of the wavelet transform in to the X values, tested in muse and simulated, works well, added test

Author: kylemath

Muse_P3example.py

@@ -6,15 +6,19 @@
 #python
 from utils import *
 data_dir = 'visual/cueing'
-subs = [101,102]
+subs = [101, 102, 103, 104, 105, 106, 108, 109, 110, 111, 112,
+        202, 203, 204, 205, 207, 208, 209, 210, 211,
+        301, 302, 303, 304, 305, 306, 307, 308, 309]
+        subs = [101, 102, 103, 104]
 nsesh = 2
 event_id = {'LeftCue': 1,'RightCue': 2}
 #Load Data
 raw = LoadMuseData(subs,nsesh,data_dir)
 #Pre-Process EEG Data
-epochs = PreProcess(raw,event_id)
+epochs = PreProcess(raw,event_id,epoch_time=(-1,2))
 #Engineer Features for Model
-feats = FeatureEngineer(epochs)
+feats = FeatureEngineer(epochs,frequency_domain=True,
+						include_phase=True)
 #Create Model
 model,_ = CreateModel(feats)
 #Train with validation, then Test

simulate.py

@@ -9,15 +9,17 @@
 # 	raw,event_id = SimulateRaw(amp1=10, amp2=5, freq=2., batch=1)
 # 	raw.save(raw_filename,overwrite=True)
 
-raw,event_id = SimulateRaw(amp1=100, amp2=50, freq=2., batch=4)
-epochs = PreProcess(raw,event_id,filter_data=False,plot_erp=True)
+raw,event_id = SimulateRaw(amp1=100, amp2=50, freq=2., batch=1)
+epochs = PreProcess(raw,event_id,filter_data=False,plot_erp=False,
+					epoch_time=(-1,2))
 
 pick = 33
 for event in event_id.keys():
 	fig = plt.imshow(epochs[event]._data[:,pick,:])
 	plt.show()
 
-feats = FeatureEngineer(epochs,model_type='CNN')
-model,_ = CreateModel(feats, units=[256,128,128,64,32,16])
+feats = FeatureEngineer(epochs,model_type='NN',
+						frequency_domain=True,include_phase=True)
+model,_ = CreateModel(feats, units=[64,32,16])
 TrainTestVal(model,feats)
 

tests.py

@@ -76,6 +76,28 @@ def test_simulate_raw(self):
 
         self.assertLess(data['acc'], 1)
 
+    def test_frequencydomain_complex(self):
+        """
+        Testing simulated data pipeline.
+        """
+        # Simulate Data
+        raw,event_id = SimulateRaw(amp1=50, amp2=60, freq=1.)
+
+        # Pre-Process EEG Data
+        epochs = PreProcess(raw,event_id)
+
+        # Engineer Features for Model
+        feats = FeatureEngineer(epochs,frequency_domain=True,
+                                include_phase=True)
+
+        # Create Model
+        model, _ = CreateModel(feats, units=[16,16])
+
+        # Train with validation, then Test
+        model, data = TrainTestVal(model,feats, 
+                    train_epochs=1,show_plots=False)
+
+        self.assertLess(data['acc'], 1)
 
 if __name__ == '__main__':
     unittest.main()

utils.py

@@ -377,8 +377,6 @@ def PreProcess(raw, event_id, plot_psd=False, filter_data=True,
                emcp_raw=False, emcp_epochs=False, epoch_decim=1, plot_electrodes=False,
                plot_erp=False):
 
-
-
   sfreq = raw.info['sfreq']
   #create new output freq for after epoch or wavelet decim
   nsfreq = sfreq/epoch_decim
@@ -459,9 +457,10 @@ def PreProcess(raw, event_id, plot_psd=False, filter_data=True,
 
 def FeatureEngineer(epochs, model_type='NN',
                     frequency_domain=False,
-                    normalization=True, electrode_median=False,
-                    wavelet_decim=1,flims=(3,30),
-                    f_bins=20,wave_cycles=3,
+                    normalization=False, electrode_median=False,
+                    wavelet_decim=1, flims=(3,30), include_phase=False,
+                    f_bins=20, wave_cycles=3, 
+                    wavelet_electrodes = [11,12,13,14,15],
                     spect_baseline=[-1,-.5],
                     test_split = 0.2, val_split = 0.2,
                     random_seed=1017, watermark = False):
@@ -492,51 +491,50 @@ def FeatureEngineer(epochs, model_type='NN',
     f_low = flims[0]
     f_high = flims[1]
     frequencies =  np.linspace(f_low, f_high, f_bins, endpoint=True)
-    eeg_chans = pick_types(epochs.info,eeg=True,eog=False)
 
-    ####
-    ## Condition0 ##
-    print('Computing Morlet Wavelets on ' + event_names[0])
-    tfr0 = tfr_morlet(epochs[event_names[0]], freqs=frequencies,
-                          n_cycles=wave_cycles, return_itc=False,
-                          picks=eeg_chans, average=False,
-                          decim=wavelet_decim)
-    tfr0 = tfr0.apply_baseline(spect_baseline,mode='mean')
-    #reshape data
-    stim_onset = np.argmax(tfr0.times>0)
-    feats.new_times = tfr0.times[stim_onset:]
-
-    #move electrodes last
-    cond0_power_out = np.moveaxis(tfr0.data[:,:,:,stim_onset:],1,3)
-    # move time second
-    cond0_power_out = np.moveaxis(cond0_power_out,1,2)
-    ####
-
-    ####
-    ## Condition1 ##
-    print('Computing Morlet Wavelets on ' + event_names[1])
-    tfr1 = tfr_morlet(epochs[event_names[1]], freqs=frequencies,
-                          n_cycles=wave_cycles, return_itc=False,
-                          picks=eeg_chans, average=False,
-                          decim=wavelet_decim)
-    tfr1 = tfr1.apply_baseline(spect_baseline,mode='mean')
-    #reshape data
-    cond1_power_out = np.moveaxis(tfr1.data[:,:,:,stim_onset:],1,3)
-    cond1_power_out = np.moveaxis(cond1_power_out,1,2) # move time second
-    ####
-
-    print('Condition one trials: ' + str(len(cond1_power_out)))
-    print(event_names[1] + ' Time Points: ' + str(len(feats.new_times)))
-    print(event_names[1] + ' Frequencies: ' + str(len(tfr1.freqs)))
-    print('Condition zero trials: ' + str(len(cond0_power_out)))
-    print(event_names[0] + ' Time Points: ' + str(len(feats.new_times)))
-    print(event_names[0] + ' Frequencies: ' + str(len(tfr0.freqs)))
+    if wavelet_electrodes == 'all':
+      wavelet_electrodes = pick_types(epochs.info,eeg=True,eog=False)
 
+    #type of output from wavelet analysis
+    if include_phase:
+      tfr_output_type = 'complex'
+    else:
+      tfr_output_type = 'power'
+
+    tfr_dict = {}
+    for event in event_names:
+      print('Computing Morlet Wavelets on ' + event)
+      tfr_temp = tfr_morlet(epochs[event], freqs=frequencies,
+                            n_cycles=wave_cycles, return_itc=False,
+                            picks=wavelet_electrodes, average=False,
+                            decim=wavelet_decim, output=tfr_output_type)
+      tfr_temp = tfr_temp.apply_baseline(spect_baseline,mode='mean')
+      stim_onset = np.argmax(tfr_temp.times>0)
+      power_out_temp = np.moveaxis(tfr_temp.data[:,:,:,stim_onset:],1,3)
+      power_out_temp = np.moveaxis(power_out_temp,1,2)
+      print(event + ' trials: ' + str(len(power_out_temp)))
+      tfr_dict[event] = power_out_temp
+
+    #reshape data (sloppy but just use the last temp tfr)
+    feats.new_times = tfr_temp.times[stim_onset:]
+
+    for event in event_names:
+      print(event + ' Time Points: ' + str(len(feats.new_times)))
+      print(event + ' Frequencies: ' + str(len(tfr_temp.freqs)))
 
     #Construct X and Y
-    X = np.append(cond0_power_out,cond1_power_out,0);
-    Y_class = np.append(np.zeros(len(cond0_power_out)),
-                        np.ones(len(cond1_power_out)),0)
+    for ievent,event in enumerate(event_names):
+      if ievent == 0:
+        X = tfr_dict[event]
+        Y_class = np.zeros(len(tfr_dict[event]))
+      else:
+        X = np.append(X,tfr_dict[event],0)
+        Y_class = np.append(Y_class,np.ones(len(tfr_dict[event]))*ievent,0)
+
+    #concatenate real and imaginary data
+    if include_phase:
+      print('Concatenating the real and imaginary components')
+      X = np.append(np.real(X),np.imag(X),2)
 
     if electrode_median:
       print('Computing Median over electrodes')