External Email - Use Caution
Of course, and thanks for looking into this. I have three pipelines the first one is short, just trying to read the epochs but gives me the epochs.times error:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import nibabel as nib
import pylab as pl
import mne
import scipy
from mne.time_frequency.tfr import cwt, morlet
import sys
sys.path.append('/N/u/mattwint/Carbonate/Desktop/net/')
import hcp
from hcp import preprocessing as preproc
subject='100307'
hcp_path='/N/dcwan/projects/hcp'
raw1=hcp.read_epochs(subject, 'task_working_memory', onset='stim', run_index=0, hcp_path=hcp_path, return_fixations_motor=False)
raw2=hcp.read_epochs(subject, 'task_working_memory', onset='stim', run_index=1, hcp_path=hcp_path, return_fixations_motor=False)
The second script is my attempt to work with the raw instance and take into account some of the things the GitHub mne-hcp page says to do, although I do not understand it all. I can run this script all they way through, but end up getting odd cardiac components I do not trust. I have tried running it with and without the reference correction and with the ecg+ and ecg- (but not both at the same time) and I have not gotten a component back that looks right to me, although I am more use to eeg data, so maybe it is just my ignorance. Note I pull the trial info, but have not gotten to the part were I epoch it yet:.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import nibabel as nib
import pylab as pl
import mne
import scipy
from mne.time_frequency.tfr import cwt, morlet
import sys
sys.path.append('/N/u/mattwint/Carbonate/Desktop/data/')
import hcp
data_type='task_working_memory'
subject='100307'
hcp_path='/N/dcwan/projects/hcp'
raw1=hcp.read_raw(subject, data_type, run_index=0, hcp_path=hcp_path)
raw2=hcp.read_raw(subject, data_type, run_index=1, hcp_path=hcp_path)
#%%########################### trial info #####################################
trial_info1=hcp.read_trial_info(subject, data_type, run_index=0, hcp_path=hcp_path)
trial_info2=hcp.read_trial_info(subject, data_type, run_index=1, hcp_path=hcp_path)
info1=len(trial_info1['stim']['codes'])
info2=len(trial_info2['stim']['codes'])
trial_info=np.zeros([info1 + info2, 40])
trial_info[0:info1,0:40]=(trial_info1['stim']['codes'])
trial_info[info1:info1+info2,0:40]=(trial_info2['stim']['codes'])
#%%
concatenated=mne.concatenate_raws([raw1,raw2])
concatenated.annotations
concatenated.annotations.delete([0,1])
concatenated.plot()
raw = concatenated.copy()
raw.load_data()
#%%###################### ??? regress out ref. ??? ############################
hcp.preprocessing.apply_ref_correction(raw)
#%%#Filter the data to 0.1-100 Hz (IIR for HP, FIR for LP, as was established previously)
raw.filter(l_freq = 0.1, h_freq = None, method = 'iir')
raw.filter(h_freq = 100, l_freq = None, method = 'fir')
#%%#View the PSD to make sure our filter took effect and orient cap to head
hcp.preprocessing.map_ch_coords_to_mne(raw)
raw.plot_psd(area_mode='range', average=False, fmax = 200)
#%%
room_noise=hcp.read_raw(subject, 'noise_empty_room', hcp_path=hcp_path)
room_info = hcp.read_info(subject, 'noise_empty_room', hcp_path=hcp_path)
#%%
empty_room_projs = mne.compute_proj_raw(room_noise, n_grad=3, n_mag=3)
raw.add_proj(empty_room_projs)
#%%
system_projs = raw.info['projs']
mags = mne.pick_types(raw.info, meg='mag')
for title, projs in [('system', system_projs),
('subject-specific', empty_room_projs[3:])]:
raw.add_proj(projs, remove_existing=True)
fig = raw.plot(proj=True, order=mags, duration=1, n_channels=2)
fig.subplots_adjust(top=0.9) # make room for title
fig.suptitle('{} projectors'.format(title), size='xx-large', weight='bold')
#%% ???(Do I need + & -)???
#The following few cells describe the process of deleting ECG projections
ecg_evoked = mne.preprocessing.create_ecg_epochs(raw,ch_name='ECG-').average()
ecg_evoked.apply_baseline((None, None))
ecg_evoked.plot_joint()
#Let's try making a projection based on our ECG data
projs, events = mne.preprocessing.compute_proj_ecg(raw, ch_name = 'ECG-', n_grad=1, n_mag=1, n_eeg=1, reject=None)
#Here we will print out the list of projectors
print(projs)
#This creates a scalp graphic with the distribution of projectors
mne.viz.plot_projs_topomap(projs, info=raw.info)
raw.add_proj(projs)
#%%#Let's do the same for vertical EOG
eog_evoked = mne.preprocessing.create_eog_epochs(raw, ch_name = 'VEOG-').average()
eog_evoked.apply_baseline((None, None))
eog_evoked.plot_joint()
#Compute projections; no_proj is set to "True" to avoid printing the whole set of projections computed above
eog_projs, _ = mne.preprocessing.compute_proj_eog(raw, n_grad=1, n_mag=1, n_eeg=1, reject=None, ch_name = 'VEOG-', no_proj=True)
#Plot the results of projections of vertical EOG
mne.viz.plot_projs_topomap(eog_projs, info=raw.info)
for title in ('Without', 'With'):
if title == 'With':
raw.add_proj(eog_projs)
fig = raw.plot()
fig.subplots_adjust(top=0.9) # make room for title
fig.suptitle('{} EOG projectors'.format(title), size='xx-large',
weight='bold')
mne.viz.plot_projs_topomap(eog_projs, info=raw.info)
raw.add_proj(eog_projs)
#%%#And horizontal EOG
eogh_evoked = mne.preprocessing.create_eog_epochs(raw, ch_name = 'HEOG-').average()
eogh_evoked.apply_baseline((None, None))
eogh_evoked.plot_joint()
#Compute projections; no_proj is set to "True" to avoid printing the whole set of projections computed above
eogh_projs, _ = mne.preprocessing.compute_proj_eog(raw, n_grad=1, n_mag=1, n_eeg=1, reject=None, ch_name = 'HEOG-', no_proj=True)
#Plot the results of projections of horizontal EOG
mne.viz.plot_projs_topomap(eogh_projs, info=raw.info)
for title in ('Without', 'With'):
if title == 'With':
raw.add_proj(eogh_projs)
fig = raw.plot()
fig.subplots_adjust(top=0.9) # make room for title
fig.suptitle('{} EOG projectors'.format(title), size='xx-large',
weight='bold')
mne.viz.plot_projs_topomap(eogh_projs, info=raw.info)
raw.add_proj(eogh_projs)
#%%
raw.apply_proj
#%%#Let's back it up one more time
raw.save('concatenated'+'filt01-100_'+'SSPapplied_tsss.fif', overwrite=True)
#%%#Filter the data to 1-40, as was established previously
raw.filter(l_freq = 1, h_freq = None, method = 'iir')
raw.filter(h_freq = 55, l_freq = None, method = 'fir')
#%%#And back it up
raw.save('concatenated'+'filt1-40_tsss.fif', overwrite=True)
#%%
raw.plot()
raw.plot_psd(fmin = 1, fmax = 40)
My third script works fine, I pull over parts of the tmegpreproc files and concatenate them back together. The only problem I run into on that one is there are a different number of good sensors between run 1 and run 2, so I have yet to make a decision on the best way to work around that problem. If you want to see that I can send it too, it just doesn?t really utilize the mne-hcp tools.