# -*- coding: utf-8 -*-
"""
---
Title: Photodiode Analysis
Date: 15-01-2018
Author: D.Cook
---
"""
import os
import json
import mne
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
Introduction¶
Testing the temporal accuracy of triggers indicating stimulus onsets¶
Event-related potentials require precise measurements of stimulus onsets. This is accomplished by sending triggers from the stimulus computer to the acquisition computer. Triggers indicate when a stimulus was shown on the EEG LAB's CRT monitor. In order to measure the temporcal accuracy of the triggers, a simple experiment was conducted with the BioSemi photodiode attached to the CRT monitor. This notebook documents the matlab code used to send the triggers and stimuli, and the python analysis of the output.
#from IPython.display import Image
Image('/Users/cookdj0128/Desktop/eegLab/photoDiodeSetup.png', retina=True)
Figure 1. The Expeirmental Setup. A photodiode is placed on the top corner at coordinate (x,y).
In MATLAB, Psychtoolbox is used to open a screen, make the screen black, and flash a white square in the top corner of the screen for $t$ seconds. Using io64.mexw64, markers are sent to BioSemi Actiview via Parallel Port to test the timing of the system. In particular, we want to know if a) the markers arrive earlier / later than the stimulus onset, b) if so, if the difference is constant / variable. The MATLAB code for the experiment is given below:
Note: The following code has been modified from the original code, which was used to analyze the current data. The only difference between the original and the below code is that a trigger code of 10 was sent to the acquistion computer prior to the Screen('Flip') command for the white box, and a code of 20 after the offset of the white box.
.m
DoTriggers = 1;
%we are not using PsychDefaultSetup
%PsychDefaultSetup(1);
% Settings for the experiment:
P.WhichScreen = 2;
P.myWidth = 1280;
P.myHeight = 1024;
P.myRate = 85;
P.SkipSyncTests = 0;
% Tell PTB about our screen
Screen('Resolution', P.WhichScreen, P.myWidth, P.myHeight, P.myRate);
% Windows OS is buggy. We can - but do not - use the following work-arounds:
%Screen('Preference','ConserveVRAM', 16384);
%Screen('Preference', 'SkipSyncTests', P.SkipSyncTests);
% Open the screen corresponding to a CRT Monitor - Samsung SyncMaster 959NF-B
[screenNumber, rect] = Screen( 'OpenWindow', P.WhichScreen);
% Set some colors:
black=BlackIndex(screenNumber);
white=WhiteIndex(screenNumber);
% Should we hide the cursor:
%HideCursor;
% Define refresh rate:
%ifi = Screen('GetFlipInterval', screenNumber);
%whenArg=ifi * 100;
% Fill the screen with black, and present it on the screen
Screen('FillRect',screenNumber, black);
Screen('Flip',screenNumber);
% Are we sending triggers, if so, do the following:
if DoTriggers; OpenTriggerPort; StartSaveBDF; SendTrigger(100, 0.001); end
WaitSecs(1);
% Begin experiment:
for i = 1:50
SendTrigger(5, 0.001);
framesSinceLastWait = Screen('WaitBlanking', screenNumber, 25);
Screen('FillRect',screenNumber,white, [0 0 200 200] );
SendTrigger(6, 0.001);
Screen('Flip',screenNumber);
SendTrigger(7, 0.001);
framesSinceLastWait = Screen('WaitBlanking', screenNumber, 25);
Screen('FillRect',screenNumber,black);
SendTrigger(8, 0.001);
Screen('Flip',screenNumber);
SendTrigger(9, 0.001);
framesSinceLastWait = Screen('WaitBlanking', screenNumber, 25);
SendTrigger(10, 0.001);
end
% Close the Parallel Port, and the screen. END.
WaitSecs(1);
if DoTriggers; SendTrigger(200, 0.001); PauseSaveBDF; CloseTriggerPort; end
Screen('Close', screenNumber);
%ShowCursorFunctions for analyzing¶
def plotTrials(test=mne.io.Raw, triggerCodes=list):
"""
pass 2 trigger codes corresponding to on/off
returns all trials in a plot
"""
test=test
triggerCodes=triggerCodes
fig = plt.figure(figsize=(14, 4))
for x in range(len(triggerCodes)):
markerCode= str(triggerCodes[x])
theEpochs = mne.Epochs(test, events, event_id=triggerCodes[x], tmin=0.0, tmax=1.0, baseline=(0.0,0), picks=[1])
#get the data from epochs
dat0 = theEpochs.get_data()
timePoints = dat0.shape[2]
#
subplotNums = [121, 122]
#plot the data for the epochs on top of one another
ax = plt.subplot(subplotNums[x]);
for z in range(theEpochs.events.shape[0]):
plt.plot(np.linspace(0, 1000.0/timePoints, num=timePoints), dat0[z, 0, 0:])
ax.set_xlim([0.0, 0.05]);
ax.set_ylim([-0.5, 2.5]);
ax.set_title("All Trials - Trigger: "+markerCode, size='x-large');
ax.set_xlabel(r'Time from trigger (s)', size='x-large');
if x == 0:
ax.set_ylabel(r'Sensor Signal', size='xx-large');
def showAverageSensor(test=mne.io.Raw, triggerCode=10):
EEG = mne.Epochs(test, events, event_id=triggerCode, tmin=-0.2, tmax=0.5, baseline=(-0.2,0))
b = EEG.average()
fig, ax = plt.subplots(1, 3, figsize=(14, 4))
b.plot(axes=ax[0], show=False);# xlim=(-10, 200),
b.plot(axes=ax[1], xlim=(-10, 20), show=False);#hline=[0.0], xlim=(-10, 50),
EEG.plot_psd(ax=ax[2], fmax=500, show=False);
for ax, title in zip(ax[:3], ['Average Sensor Signal\n Trig: '+str(triggerCode), 'Average Sensor Signal: '+str(triggerCode), 'PSD for Sensor']):
ax.set_title(title)
plt.show()
Read the data¶
After collecting the data from the experiment, MNE is used to analyze the results.
#test = mne.io.read_raw_edf("/Volumes/MYSTUFF2016/EEGLabTriggers/jan09_photoDiode10.bdf")
test = mne.io.read_raw_edf("/Users/cookdj0128/Desktop/eegLab/testNov23_photoDiode3.bdf")
events = mne.find_events(test)
Show all trials¶
plotTrials(test, triggerCodes=[10,20])
Show Average Sensor Signal from .m line before Screen flip¶
showAverageSensor(test, 10)
Show Average Sensor Signal from .m line after Screen flip¶
showAverageSensor(test, 20)
Conclusion¶
On average, the stimulus appears after about 1 refresh rate on the CRT monitor with a resolution of 1280 x 2014 @ 85hz. Researchers are advised to adjust their analysis pipelines to reflect the refresh rate constant. Alternatively, researchers are advised to use the photodiode with a small flashing white box appearing at the same time as their stimuli to measure the precise onset of the stimulus.*
theEpochs = mne.Epochs(test, events, event_id=10, tmin=0.0, tmax=1.0, baseline=(0.0,0), picks=[1])
b = theEpochs.average()
fig, ax = plt.subplots(1, 1, figsize=(14, 4))
b.plot(axes=ax, show=False, xlim=(-10, 20));
plt.show()
theEpochs = mne.Epochs(test, events, event_id=20, tmin=0.0, tmax=1.0, baseline=(0.0,0), picks=[1])
b = theEpochs.average()
fig, ax = plt.subplots(1, 1, figsize=(14, 4))
b.plot(axes=ax, show=False, xlim=(-10, 20));
plt.show()
Figure 3. The average onset (top) and offset (bottom) of the stimulus comes after about 1 refresh rate on the CRT monitor @ 85hz or 11.76ms.