Find Photodiode Events

In this example, we use pd-parser to find photodiode events and align them to behavior. Then, we save the data to BIDS format.

# Authors: Alex Rockhill <aprockhill@mailbox.org>
#
# License: BSD (3-clause)

Simulate data and use it to make a raw object

We’ll make an mne.io.Raw object so that we can save out some random data with a photodiode event channel in it in fif format (a commonly used electrophysiology data format).

import os.path as op
import numpy as np

import mne
from mne.utils import _TempDir

import pd_parser
from pd_parser.parse_pd import _read_raw, _to_tsv

out_dir = _TempDir()
print(f'After running this example, you can find the data here: {out_dir}')

# simulate photodiode data
n_events = 300
prop_corrupted = 0.01
raw, beh, events, corrupted_indices = \
    pd_parser.simulate_pd_data(n_events=n_events,
                               prop_corrupted=prop_corrupted)

# make fake electrophysiology data
info = mne.create_info(['ch1', 'ch2', 'ch3'], raw.info['sfreq'],
                       ['seeg'] * 3)
raw2 = mne.io.RawArray(np.random.random((3, raw.times.size)) * 1e-6, info)
raw.add_channels([raw2])
raw.info['line_freq'] = 60  # for bids

fname = op.join(out_dir, 'sub-1_task-mytask_raw.fif')
raw.save(fname)

# roundtrip so that raw is properly loaded from disk and has a filename
raw = _read_raw(fname)

After running this example, you can find the data here: /tmp/tmp_mne_tempdir_2_ee581j
Creating RawArray with float64 data, n_channels=1, n_times=2042106
    Range : 0 ... 2042105 =      0.000 ...  2042.105 secs
Ready.
Creating RawArray with float64 data, n_channels=3, n_times=2042106
    Range : 0 ... 2042105 =      0.000 ...  2042.105 secs
Ready.
Writing /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif
Closing /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif
[done]
Reading in /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif
Opening raw data file /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif...
    Range : 0 ... 2042105 =      0.000 ...  2042.105 secs
Ready.

Make behavior data

We’ll make a dictionary with lists for the events that are time-stamped when the photodiode was turned on and other events relative to those events. We’ll add some noise to the time-stamps so that we can see how behavior might look in an experimental setting. Let’s make a task where there is a fixation stimulus, then a go cue, and a then response as an example.

np.random.seed(12)
# add some noise to make it harder to align, use just over
# the exclusion of 0.03 to make some events excluded
offsets = np.random.random(n_events) * 0.035 - 0.0125
# in this example, the fixation would always be 700 ms
# after which point a cue would appear which is the "go time"
go_time = np.repeat(0.7, n_events)
# let's make the response time between 0.5 and 1.5 seconds uniform random
response_time = list(go_time + np.random.random(n_events) + 1.5)
for i in [10, 129, 232, 288]:
    response_time[i] = 'n/a'  # make some no responses


# put in dictionary to be converted to tsv file
beh['fix_onset_time'] = beh['time'] + offsets
beh['go_time'] = go_time
beh['response_time'] = response_time
behf = op.join(out_dir, 'sub-1_task-mytask_beh.tsv')
# save behavior file out
_to_tsv(behf, beh)

Use the interactive graphical user interface (GUI) to find parameters

On the documentation webpage, this is example is not interactive, but if you download it as a jupyter notebook and run it or copy the code into a console running python (ipython recommended), you can see how to interact with the photodiode data to pick reasonable parameters by following the instructions.

pd_parser.find_pd_params(raw, pd_ch_names=['pd'])

Reading 0 ... 2042105  =      0.000 ...  2042.105 secs...

Find the photodiode events relative to the behavioral timing of interest

This function will use the default parameters or the parameters you found from pd_parser.find_pd_parameters() to find and align the photodiode events, excluding events that were off because the computer hung up on computation, for instance. That data is saved in the same folder as the raw file (in this case, a temperary directory generated by _TempDir()). The data can be used directly, or it can be accessed via pd_parser.pd_parser_save_to_bids() to store it in the brain imagine data structure (BIDS) standardized format before using it.

pd_parser.parse_pd(raw, beh=beh, pd_ch_names=['pd'], max_len=1.5)

• 
• 

Reading 0 ... 2042105  =      0.000 ...  2042.105 secs...
Finding photodiode events

  0%|          | 0/793 [00:00<?, ?it/s]
 77%|#######6  | 607/793 [00:00<00:00, 6065.00it/s]
100%|##########| 793/793 [00:00<00:00, 6431.41it/s]
299 up-deflection photodiode candidate events found
Checking best alignments

  0%|          | 0/300 [00:00<?, ?it/s]
  1%|1         | 3/300 [00:00<00:10, 27.80it/s]
  2%|2         | 6/300 [00:00<00:10, 27.93it/s]
  3%|3         | 9/300 [00:00<00:10, 27.74it/s]
  4%|4         | 12/300 [00:00<00:10, 27.69it/s]
  5%|5         | 15/300 [00:00<00:10, 27.66it/s]
  6%|6         | 18/300 [00:00<00:10, 27.67it/s]
  7%|7         | 21/300 [00:00<00:10, 27.41it/s]
  8%|8         | 24/300 [00:00<00:10, 27.42it/s]
  9%|9         | 27/300 [00:00<00:10, 27.28it/s]
 10%|#         | 30/300 [00:01<00:09, 27.38it/s]
 11%|#1        | 33/300 [00:01<00:09, 27.20it/s]
 12%|#2        | 36/300 [00:01<00:09, 27.23it/s]
 13%|#3        | 39/300 [00:01<00:09, 27.16it/s]
 14%|#4        | 42/300 [00:01<00:09, 27.01it/s]
 15%|#5        | 45/300 [00:01<00:09, 26.43it/s]
 16%|#6        | 48/300 [00:01<00:09, 26.71it/s]
 17%|#7        | 51/300 [00:01<00:09, 26.93it/s]
 18%|#8        | 54/300 [00:01<00:09, 26.79it/s]
 19%|#9        | 57/300 [00:02<00:09, 26.80it/s]
 20%|##        | 60/300 [00:02<00:09, 25.69it/s]
 21%|##1       | 63/300 [00:02<00:09, 25.77it/s]
 22%|##2       | 66/300 [00:02<00:09, 25.69it/s]
 23%|##3       | 69/300 [00:02<00:08, 25.81it/s]
 24%|##4       | 72/300 [00:02<00:08, 26.30it/s]
 25%|##5       | 75/300 [00:02<00:08, 26.70it/s]
 26%|##6       | 78/300 [00:02<00:08, 26.81it/s]
 27%|##7       | 81/300 [00:03<00:08, 27.01it/s]
 28%|##8       | 84/300 [00:03<00:08, 26.76it/s]
 29%|##9       | 87/300 [00:03<00:08, 25.96it/s]
 30%|###       | 90/300 [00:03<00:08, 26.18it/s]
 31%|###1      | 93/300 [00:03<00:07, 26.02it/s]
 32%|###2      | 96/300 [00:03<00:08, 25.37it/s]
 33%|###3      | 99/300 [00:03<00:08, 24.70it/s]
 34%|###4      | 102/300 [00:03<00:07, 25.14it/s]
 35%|###5      | 105/300 [00:03<00:07, 25.58it/s]
 36%|###6      | 108/300 [00:04<00:07, 25.72it/s]
 37%|###7      | 111/300 [00:04<00:07, 25.98it/s]
 38%|###8      | 114/300 [00:04<00:07, 26.13it/s]
 39%|###9      | 117/300 [00:04<00:06, 26.20it/s]
 40%|####      | 120/300 [00:04<00:06, 26.01it/s]
 41%|####1     | 123/300 [00:04<00:06, 25.93it/s]
 42%|####2     | 126/300 [00:04<00:06, 26.03it/s]
 43%|####3     | 129/300 [00:04<00:06, 26.13it/s]
 44%|####4     | 132/300 [00:04<00:06, 26.37it/s]
 45%|####5     | 135/300 [00:05<00:06, 26.37it/s]
 46%|####6     | 138/300 [00:05<00:06, 26.14it/s]
 47%|####6     | 141/300 [00:05<00:06, 26.24it/s]
 48%|####8     | 144/300 [00:05<00:06, 25.94it/s]
 49%|####9     | 147/300 [00:05<00:05, 25.64it/s]
 50%|#####     | 150/300 [00:05<00:05, 25.70it/s]
 51%|#####1    | 153/300 [00:05<00:05, 25.72it/s]
 52%|#####2    | 156/300 [00:05<00:05, 25.68it/s]
 53%|#####3    | 159/300 [00:06<00:05, 24.69it/s]
 54%|#####4    | 162/300 [00:06<00:05, 24.91it/s]
 55%|#####5    | 165/300 [00:06<00:05, 25.08it/s]
 56%|#####6    | 168/300 [00:06<00:05, 25.05it/s]
 57%|#####6    | 171/300 [00:06<00:05, 25.22it/s]
 58%|#####8    | 174/300 [00:06<00:05, 24.44it/s]
 59%|#####8    | 177/300 [00:06<00:05, 24.33it/s]
 60%|######    | 180/300 [00:06<00:04, 24.80it/s]
 61%|######1   | 183/300 [00:07<00:04, 25.21it/s]
 62%|######2   | 186/300 [00:07<00:04, 25.39it/s]
 63%|######3   | 189/300 [00:07<00:04, 25.35it/s]
 64%|######4   | 192/300 [00:07<00:04, 25.24it/s]
 65%|######5   | 195/300 [00:07<00:04, 25.39it/s]
 66%|######6   | 198/300 [00:07<00:04, 25.37it/s]
 67%|######7   | 201/300 [00:07<00:04, 24.45it/s]
 68%|######8   | 204/300 [00:07<00:03, 24.67it/s]
 69%|######9   | 207/300 [00:07<00:03, 24.88it/s]
 70%|#######   | 210/300 [00:08<00:03, 24.96it/s]
 71%|#######1  | 213/300 [00:08<00:03, 25.07it/s]
 72%|#######2  | 216/300 [00:08<00:03, 25.23it/s]
 73%|#######3  | 219/300 [00:08<00:03, 25.39it/s]
 74%|#######4  | 222/300 [00:08<00:03, 25.04it/s]
 75%|#######5  | 225/300 [00:08<00:02, 25.12it/s]
 76%|#######6  | 228/300 [00:08<00:02, 25.13it/s]
 77%|#######7  | 231/300 [00:08<00:02, 25.07it/s]
 78%|#######8  | 234/300 [00:09<00:02, 24.94it/s]
 79%|#######9  | 237/300 [00:09<00:02, 25.12it/s]
 80%|########  | 240/300 [00:09<00:02, 24.95it/s]
 81%|########1 | 243/300 [00:09<00:02, 24.89it/s]
 82%|########2 | 246/300 [00:09<00:02, 24.85it/s]
 83%|########2 | 249/300 [00:09<00:02, 24.80it/s]
 84%|########4 | 252/300 [00:09<00:02, 23.54it/s]
 85%|########5 | 255/300 [00:09<00:01, 23.40it/s]
 86%|########6 | 258/300 [00:10<00:01, 23.32it/s]
 87%|########7 | 261/300 [00:10<00:01, 21.99it/s]
 88%|########8 | 264/300 [00:10<00:01, 20.20it/s]
 89%|########9 | 267/300 [00:10<00:01, 19.39it/s]
 90%|######### | 270/300 [00:10<00:01, 19.36it/s]
 91%|#########1| 273/300 [00:10<00:01, 19.89it/s]
 92%|#########2| 276/300 [00:10<00:01, 20.38it/s]
 93%|#########3| 279/300 [00:11<00:01, 20.84it/s]
 94%|#########3| 282/300 [00:11<00:00, 21.23it/s]
 95%|#########5| 285/300 [00:11<00:00, 20.93it/s]
 96%|#########6| 288/300 [00:11<00:00, 21.11it/s]
 97%|#########7| 291/300 [00:11<00:00, 21.48it/s]
 98%|#########8| 294/300 [00:11<00:00, 20.78it/s]
 99%|#########9| 297/300 [00:11<00:00, 20.74it/s]
100%|##########| 300/300 [00:12<00:00, 21.22it/s]
100%|##########| 300/300 [00:12<00:00, 24.76it/s]
Best alignment is with the first behavioral event shifted 0.00 s relative to the first synchronization event and has errors: min -31.60 ms, q1 -9.87 ms, med 0.17 ms, q3 9.68 ms, max 34.10 ms, 3 missed events
Excluding events that have zero close synchronization events or more than one synchronization event within `max_len` time
1 off by 34 ms
57 none found
76 off by -30 ms
96 none found
97 off by 32 ms
98 off by -31 ms
172 none found

(<Annotations | 293 segments: Fixation (293)>, [12270, 'n/a', 24994, 32081, 38711, 45439, 51458, 58189, 65602, 72878, 79973, 86136, 93477, 100763, 107011, 113960, 119991, 126166, 132641, 138878, 146016, 153243, 159760, 166238, 172405, 178531, 185600, 192499, 198583, 205303, 211906, 219178, 226255, 233158, 239987, 247411, 254891, 261398, 267758, 274953, 281049, 287596, 293701, 300180, 306286, 312721, 319906, 327264, 334453, 341296, 348220, 354762, 361015, 367669, 374768, 380862, 386893, 'n/a', 400499, 407551, 414653, 422052, 428652, 435190, 442400, 449547, 456526, 463742, 470705, 478141, 484642, 491749, 499174, 505676, 512593, 519142, 'n/a', 531947, 537976, 545115, 551794, 558728, 565838, 572466, 579017, 585271, 592462, 599712, 606813, 614126, 621098, 628142, 634432, 641583, 648084, 654741, 'n/a', 'n/a', 'n/a', 681922, 688399, 694726, 701993, 708446, 715103, 721240, 727691, 733820, 740336, 747400, 754861, 761799, 768203, 775188, 782128, 789366, 796167, 803534, 810157, 816586, 823365, 830749, 837835, 844576, 851322, 858295, 864967, 871793, 879250, 886649, 893910, 900101, 906936, 913486, 920854, 927325, 934721, 941079, 947710, 955150, 962043, 969441, 976249, 982591, 989660, 996442, 1003200, 1009405, 1015622, 1023016, 1030450, 1037899, 1044436, 1050505, 1057336, 1064440, 1071412, 1077491, 1083840, 1090899, 1096982, 1104252, 1111736, 1118108, 1124189, 1131659, 1138810, 1145491, 1152208, 1159581, 1166115, 1172195, 'n/a', 1185585, 1192607, 1200029, 1206570, 1212940, 1219019, 1226266, 1232384, 1239812, 1246871, 1254238, 1261460, 1267746, 1274010, 1281452, 1287889, 1295075, 1301331, 1307780, 1314790, 1321163, 1328291, 1335610, 1342091, 1348539, 1354581, 1361215, 1367313, 1374175, 1381207, 1388059, 1394965, 1401864, 1409218, 1416015, 1423307, 1429835, 1436046, 1442379, 1449569, 1455914, 1463073, 1470427, 1477244, 1483459, 1490394, 1496963, 1504326, 1510661, 1517012, 1523736, 1530476, 1536758, 1543320, 1550322, 1557028, 1563118, 1569741, 1576365, 1583036, 1589072, 1595483, 1602909, 1610234, 1617645, 1624783, 1631434, 1638840, 1645291, 1652181, 1659062, 1665723, 1672099, 1678595, 1685559, 1691985, 1698504, 1705035, 1711522, 1718765, 1725584, 1731612, 1738724, 1745477, 1752645, 1759387, 1765853, 1773112, 1779372, 1785789, 1792515, 1799805, 1805978, 1812486, 1819284, 1826302, 1832460, 1839593, 1846428, 1853204, 1859775, 1866247, 1873583, 1881074, 1887337, 1893654, 1900655, 1907802, 1915036, 1921490, 1927804, 1933913, 1940477, 1947321, 1954630, 1960680, 1966814, 1973794, 1980801, 1986987, 1993132, 1999174, 2006147, 2012567, 2018928, 2025925, 2032107])

Add events relative to the photodiode events

The photodiode is usually sychronized to one event (e.g. the fixation so that if the deflections caused by the photodiode are large enough to influence other channels through amplifier interactions it doesn’t cause issues with the analysis) so often the events of interest are relative to the photodiode event. In the task a timer can be started at the photodiode event and checked each time a subsequent event occurs. These events should then be recorded in tsv file, which can be passed to pd-parser in order to add the events. Note: if more than one photodiode event is used, the parser can be used for each event separately using the keyword add_event=True.

pd_parser.add_relative_events(
    raw, beh,
    relative_event_keys=['go_time', 'response_time'],
    relative_event_names=['Go Cue', 'Response'])

Overwriting existing file.

<Annotations | 875 segments: Fixation (293), Go Cue (293), Response (289)>

Save data to BIDS format

This saves our data to BIDS format so that it’s ready to be analyzed in a reproducible way; BIDS requires all the files the BIDS community has deemed necessary for analysis, so you should have everything you need to continue on with an analysis at this point. See https://bids.neuroimaging.io/ and https://mne.tools/mne-bids/ for more information about BIDS.

pd_parser.save_to_bids(op.join(out_dir, 'bids_dir'), fname,
                       sub='1', task='mytask')

Reading in /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif
Opening raw data file /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif...
    Range : 0 ... 2042105 =      0.000 ...  2042.105 secs
Ready.
Used Annotations descriptions: ['Fixation', 'Go Cue', 'Response']
Opening raw data file /tmp/tmp_mne_tempdir_2_ee581j/sub-1_task-mytask_raw.fif...
    Range : 0 ... 2042105 =      0.000 ...  2042.105 secs
Ready.
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/README'...
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/participants.tsv'...
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/participants.json'...
Used Annotations descriptions: ['Fixation', 'Go Cue', 'Response']
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/sub-1/ieeg/sub-1_task-mytask_events.tsv'...
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/dataset_description.json'...
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/sub-1/ieeg/sub-1_task-mytask_ieeg.json'...
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/sub-1/ieeg/sub-1_task-mytask_channels.tsv'...
/home/alexrockhill/software/anaconda3/envs/mne/lib/python3.10/site-packages/mne_bids/write.py:1848: RuntimeWarning: Converting data files to BrainVision format
  warn('Converting data files to BrainVision format')
Writing '/tmp/tmp_mne_tempdir_2_ee581j/bids_dir/sub-1/sub-1_scans.tsv'...
Wrote /tmp/tmp_mne_tempdir_2_ee581j/bids_dir/sub-1/sub-1_scans.tsv entry with ieeg/sub-1_task-mytask_ieeg.vhdr.