Resting State Processing Pipeline
Author: Hadi Zaatiti hadi.zaatiti@nyu.edu
In this script we will be processing data generated from the KIT-MEG for a resting state experiment.
The experiment has been conducted using the code in Link Text involving a 5 minutes eyes closed followed by a 5 minutes eyes open of the participant.
The data can be accessed at Data Storage we will be using the resting_state/sub-01 BIDS dataset
After downloading the data, be familiar with the different files in the meg-kit folder:
The .con are the raw files produced by the KIT machine, has the MEG time series for each channel
The .con that has analysis in their name, have already applied a filter to account for the noisy magnetic field in the MSR, the latter is measured with the magnetometers of the KIT, this file contains the filtered MEG time series for each channel
The .mrk has the markers position from KIT
The .fif is produced by applying KIT2FIFF command provided by the MNE environment from the filtered .con file and the .mrk files (Check the KIT2FIF tutorial in this documentation)
We obtained two .fif files, one for eyes closed and one for eyes open
Import mne and set visualisations to show in the notebook
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import mne
Load your .fif file for eyes closed
[2]:
# Load your FIFF file
raw = mne.io.read_raw_fif(r"C:\Users\hz3752\Box\MEG\Data\resting-state\sub-01\meg-kit\sub-01_01-eyes-closed-raw.fif", verbose=False)
Display the info structure from the data
[3]:
print(raw.info)
<Info | 13 non-empty values
bads: []
ch_names: MEG 001, MEG 002, MEG 003, MEG 004, MEG 005, MEG 006, MEG 007, ...
chs: 207 Magnetometers, 17 Reference Magnetometers, 32 misc, 1 Stimulus
custom_ref_applied: False
description: New York University Abu Dhabi/224-channel MEG System (442) ...
dev_head_t: MEG device -> head transform
dig: 3459 items (3 Cardinal, 5 HPI, 3451 Extra)
file_id: 4 items (dict)
highpass: 0.0 Hz
kit_system_id: 442 (New York University Abu Dhabi, 2014-)
lowpass: 500.0 Hz
meas_date: 2024-04-19 09:05:59 UTC
meas_id: 4 items (dict)
nchan: 257
projs: []
sfreq: 1000.0 Hz
>
Let us plot the sensor layout of the system used in the data acquisition (KIT-MEG)
[4]:
# For a 2D topographic plot of the sensor locations
raw.plot_sensors(kind='topomap', show_names=True);
[5]:
# For a 3D plot, you can also do:
fig = raw.plot_sensors(kind='3d', show_names=True)
Let us visualise the data and plot the first 5 seconds of the data interactively)
[6]:
%matplotlib inline
# Plot the first 5 seconds of the data
raw.plot(start=0, duration=5)
Using qt as 2D backend.
[6]:
<mne_qt_browser._pg_figure.MNEQtBrowser at 0x21f2a46fac0>
You should see the following interactive window
From the above interactive window, you can mark channels as bad by visual inspection
Notice that in MNE, the channels have different types
MEG, for our system a .fif will show
MEG 001toMEG 208are of typemagMISC, for our system this will show from
MISC 001toMISC 032are of typemiscSTIM, for our system this will show as
STI 014are of typestim
Let us print the channel types from our raw data.
[7]:
print(raw.info.get_channel_types())
['mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'ref_meg', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'mag', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'ref_meg', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'misc', 'stim']
Let us print all the channel names
[8]:
print(raw.ch_names)
['MEG 001', 'MEG 002', 'MEG 003', 'MEG 004', 'MEG 005', 'MEG 006', 'MEG 007', 'MEG 008', 'MEG 009', 'MEG 010', 'MEG 011', 'MEG 012', 'MEG 013', 'MEG 014', 'MEG 015', 'MEG 016', 'MEG 017', 'MEG 018', 'MEG 019', 'MEG 020', 'MEG 021', 'MEG 022', 'MEG 023', 'MEG 024', 'MEG 025', 'MEG 026', 'MEG 027', 'MEG 028', 'MEG 029', 'MEG 030', 'MEG 031', 'MEG 032', 'MEG 033', 'MEG 034', 'MEG 035', 'MEG 036', 'MEG 037', 'MEG 038', 'MEG 039', 'MEG 040', 'MEG 041', 'MEG 042', 'MEG 043', 'MEG 044', 'MEG 045', 'MEG 046', 'MEG 047', 'MEG 048', 'MEG 049', 'MEG 050', 'MEG 051', 'MEG 052', 'MEG 053', 'MEG 054', 'MEG 055', 'MEG 056', 'MEG 057', 'MEG 058', 'MEG 059', 'MEG 060', 'MEG 061', 'MEG 062', 'MEG 063', 'MEG 064', 'MEG 065', 'MEG 066', 'MEG 067', 'MEG 068', 'MEG 069', 'MEG 070', 'MEG 071', 'MEG 072', 'MEG 073', 'MEG 074', 'MEG 075', 'MEG 076', 'MEG 077', 'MEG 078', 'MEG 079', 'MEG 080', 'MEG 081', 'MEG 082', 'MEG 083', 'MEG 084', 'MEG 085', 'MEG 086', 'MEG 087', 'MEG 088', 'MEG 089', 'MEG 090', 'MEG 091', 'MEG 092', 'MEG 093', 'MEG 094', 'MEG 095', 'MEG 096', 'MEG 097', 'MEG 098', 'MEG 099', 'MEG 100', 'MEG 101', 'MEG 102', 'MEG 103', 'MEG 104', 'MEG 105', 'MEG 106', 'MEG 107', 'MEG 108', 'MEG 109', 'MEG 110', 'MEG 111', 'MEG 112', 'MEG 113', 'MEG 114', 'MEG 115', 'MEG 116', 'MEG 117', 'MEG 118', 'MEG 119', 'MEG 120', 'MEG 121', 'MEG 122', 'MEG 123', 'MEG 124', 'MEG 125', 'MEG 126', 'MEG 127', 'MEG 128', 'MEG 129', 'MEG 130', 'MEG 131', 'MEG 132', 'MEG 133', 'MEG 134', 'MEG 135', 'MEG 136', 'MEG 137', 'MEG 138', 'MEG 139', 'MEG 140', 'MEG 141', 'MEG 142', 'MEG 143', 'MEG 144', 'MEG 145', 'MEG 146', 'MEG 147', 'MEG 148', 'MEG 149', 'MEG 150', 'MEG 151', 'MEG 152', 'MEG 153', 'MEG 154', 'MEG 155', 'MEG 156', 'MEG 157', 'MEG 158', 'MEG 159', 'MEG 160', 'MEG 161', 'MEG 162', 'MEG 163', 'MEG 164', 'MEG 165', 'MEG 166', 'MEG 167', 'MEG 168', 'MEG 169', 'MEG 170', 'MEG 171', 'MEG 172', 'MEG 173', 'MEG 174', 'MEG 175', 'MEG 176', 'MEG 177', 'MEG 178', 'MEG 179', 'MEG 180', 'MEG 181', 'MEG 182', 'MEG 183', 'MEG 184', 'MEG 185', 'MEG 186', 'MEG 187', 'MEG 188', 'MEG 189', 'MEG 190', 'MEG 191', 'MEG 192', 'MEG 193', 'MEG 194', 'MEG 195', 'MEG 196', 'MEG 197', 'MEG 198', 'MEG 199', 'MEG 200', 'MEG 201', 'MEG 202', 'MEG 203', 'MEG 204', 'MEG 205', 'MEG 206', 'MEG 207', 'MEG 208', 'MEG 209', 'MEG 210', 'MEG 211', 'MEG 212', 'MEG 213', 'MEG 214', 'MEG 215', 'MEG 216', 'MEG 217', 'MEG 218', 'MEG 219', 'MEG 220', 'MEG 221', 'MEG 222', 'MEG 223', 'MEG 224', 'MISC 001', 'MISC 002', 'MISC 003', 'MISC 004', 'MISC 005', 'MISC 006', 'MISC 007', 'MISC 008', 'MISC 009', 'MISC 010', 'MISC 011', 'MISC 012', 'MISC 013', 'MISC 014', 'MISC 015', 'MISC 016', 'MISC 017', 'MISC 018', 'MISC 019', 'MISC 020', 'MISC 021', 'MISC 022', 'MISC 023', 'MISC 024', 'MISC 025', 'MISC 026', 'MISC 027', 'MISC 028', 'MISC 029', 'MISC 030', 'MISC 031', 'MISC 032', 'STI 014']
Let us plot the MISC 001 channel that contains the trigger, it is a good practice to copy the raw data before picking a specific channel since the picking operation changes
[9]:
%matplotlib inline
channel_name = 'MISC 001'
raw_picked = raw.copy().pick_channels([channel_name])
scalings = {'misc':0.1}
raw_picked.plot(scalings = scalings, duration=315, start=0, n_channels=1)
NOTE: pick_channels() is a legacy function. New code should use inst.pick(...).
[9]:
<mne_qt_browser._pg_figure.MNEQtBrowser at 0x21f2bdf3ac0>
You should see the trigger channel going from 0 to 1 over a period of 5 minutes, corresponding to the eyes closed period as in the following image. The AU unit on the Y-axis correspond to arbitrary unit.
Note the times of beginning and end of the trigger as respectively 11 seconds and 310 seconds, let’s do FFT on the data for this time
[10]:
croped_data = raw.copy()
[11]:
croped_data.crop(11, 310)
[11]:
General
| Measurement date | April 19, 2024 09:05:59 GMT |
|---|---|
| Experimenter | Unknown |
| Participant | Unknown |
Channels
| Digitized points | 3459 points |
|---|---|
| Good channels | 207 Magnetometers, 17 Reference Magnetometers, 32 misc, 1 Stimulus |
| Bad channels | None |
| EOG channels | Not available |
| ECG channels | Not available |
Data
| Sampling frequency | 1000.00 Hz |
|---|---|
| Highpass | 0.00 Hz |
| Lowpass | 500.00 Hz |
| Filenames | sub-01_01-eyes-closed-raw.fif |
| Duration | 00:04:59 (HH:MM:SS) |
[12]:
%matplotlib inline
croped_data.plot(start=0, duration =5)
[12]:
<mne_qt_browser._pg_figure.MNEQtBrowser at 0x21f2cf97640>
[13]:
data, times = raw[:]
wsize=256
fourier_transform_data = mne.time_frequency.stft(data, wsize=wsize)
Number of frequencies: 129
Number of time steps: 2516
[14]:
%matplotlib inline
# Select the channel index you're interested in
channel_index = 1
import numpy as np
# Compute magnitude of STFT results
magnitude = np.abs(fourier_transform_data[channel_index])
# Prepare frequencies and time vectors for plotting
sfreq = raw.info['sfreq']
print('Sampling frequency', sfreq)
frequencies = np.linspace(0, sfreq / 2, magnitude.shape[0], endpoint=True)
t_secs = np.arange(magnitude.shape[1]) * (wsize / sfreq)
# Plotting
plt.figure()
plt.pcolormesh(t_secs, frequencies, magnitude, shading='auto')
plt.title(f'STFT Magnitude - Channel {raw.info["ch_names"][channel_index]}')
plt.ylabel('Frequency (Hz)')
plt.ylim(0, 20)
plt.xlabel('Time (sec)')
plt.colorbar(label='Magnitude')
plt.show()
Sampling frequency 1000.0
[15]:
import numpy as np
import matplotlib.pyplot as plt
# Assuming 'stft_result' from previous STFT calculation and 'raw' object are available
# Parameters
channel_index = 0
sfreq = raw.info['sfreq']
wsize = 256 # Window size used in STFT
hop_size = int(wsize / 2) # 50% overlap
# Compute magnitude of STFT results
magnitude = np.abs(fourier_transform_data[channel_index])
# Prepare frequencies and time vectors for plotting
frequencies = np.linspace(0, sfreq / 2, magnitude.shape[0], endpoint=True)
# Calculate correct start times for each window and then add half the window size to center it
start_times = np.arange(0, magnitude.shape[1] * hop_size, hop_size)
t_secs = (start_times + wsize / 2) / sfreq
# Limit frequency data to 100 Hz
freq_limit = 20
freq_indices = frequencies <= freq_limit
frequencies = frequencies[freq_indices]
magnitude = magnitude[freq_indices, :]
# Plotting
plt.figure()
plt.pcolormesh(t_secs, frequencies, magnitude, shading='auto')
plt.title(f'STFT Magnitude - Channel {raw.info["ch_names"][channel_index]}')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (sec)')
plt.ylim(0, freq_limit) # Set y-axis limit to 100 Hz
plt.colorbar(label='Magnitude')
plt.show()
[ ]: