If you use the EP Toolkit for ERP analysis, you can make great topographical plots of the scalp distribution of ERPs for single points in time. Problem is that reviewers like to see plots for averages over windows of time (in, say, 100ms increments). If you want to make such plots (like the one below), you can use the following code. You need to save a grand average .ept file first, but everything else should be pretty self explanatory. Copy the code below into the Matlab editor and save it as something like ‘plotAvgTopos.m’.
function plotAvgTopos
% Function to plot average topographical maps from ERP PCA toolbox files.
%
% Depends on the EEGlab toolbox for core functionality. If you have the ERP
% PCA toolbox installed, then you already have the EEGlab toolbox as well.
%
% Requires that you have saved out a grand average .ept file using the ERP
% PCA toolbox.
%
% The resulting matlab figure can be manipulated in the figure window or
% can be saved out as an image file to be edited using an image editing
% package (e.g., photoshop).
%
% Brock Kirwan - kirwan@byu.edu
% 5/13/2013
%%Preference variables
azimuth = 80; %azimuth (rotation around the IS axis) for 3D plot
elevation = 45; %elevation (rotation around the LR axis) for 3D plot
%%Locate and load the grand average ept file
%ask the user to indicate the file
[filename, pathname] = uigetfile('*.ept','Select Grand Average EPT File');
%then load it. This loads in a structure called EPdata.
load([pathname filename],'-mat');
%%Ask user for epochs to average
windowTimes = inputdlg({'Start avg time (in ms relative to stim onset)', 'Stop avg time', 'Number of windows'}, 'Epochs');
%error checking
if str2double(windowTimes{1}) >= str2double(windowTimes{2})
error('Start time must be less than end time.');
end
%These are the windows (in ms relative to stim onset) that you want to plot
% should look something like: [300, 400, 500, 600]
windows = linspace(str2double(windowTimes{1}), ...
str2double(windowTimes{2}), ...
str2double(windowTimes{3})+1);
%%Ask user for which cells to plot
cells= listdlg('PromptString','Select conditions to plot:',...
'SelectionMode','multiple',...
'ListString',EPdata.cellNames);
%%Read in some info from the EPdata structure
%Figure out the sampling period. Usually 4ms.
period = 1/EPdata.Fs*1000;
%get the sample numbers based on the sampling period.
windInd = windows/period;
%add the baseline samples. Usually 50.
windInd = windInd + EPdata.baseline;
%Count up the number of windows and conditions to plot.
nWindows = length(windows)-1;
nCells = length(cells);
for cel = 1:nCells
cellNames{cel} = EPdata.cellNames{cells(cel)};
end
%%Ask user for and read electrode location file
if exist('/Applications/EP_Toolkit/electrodes/GSN-Hydrocel-129.ced','file');
eloc = readlocs('/Applications/EP_Toolkit/electrodes/GSN-Hydrocel-129.ced');
else
[filename, pathname] = uigetfile('*.ced','Select Electrode Location File');
eloc = readlocs([pathname filename]);
end
%get rid of FID electrodes
eloc(1:3) = [];
%%ask the user if there are any comparisons
ButtonName = questdlg('Any subtractions?', ...
'Question', ...
'Yes', 'No', 'No');
if strcmp(ButtonName,'Yes')
posCell = listdlg('PromptString', 'Select positive (+) condition:', ...
'SelectionMode', 'single', ...
'ListString', cellNames);
negCell = listdlg('PromptString','Select negative (-) condition:', ...
'SelectionMode', 'single', ...
'ListString', cellNames);
cellNames{nCells + 1} = [cellNames{posCell} '-' cellNames{negCell}];
end
%%Loop over windows and conditions to get averaged data
%place-holder variable
data = zeros(size(EPdata.data,1),nWindows,nCells);
for win = 1:nWindows
for cel = 1:nCells
data(:,win,cel) = mean(EPdata.data(:,windInd(win):windInd(win+1),cells(cel)),2);
end
end
if strcmp(ButtonName,'Yes')
data = cat(3,data,data(:,:,posCell) - data(:,:,negCell));
nCells = nCells + 1;
end
%%set the min/max scale for the color map
dataMin = min(min(min(data)));
dataMax = max(max(max(data)));
prompt={'Enter the minimum value:','Enter the maximum value:'};
name='Values for color scale';
numlines=1;
defaultanswer={num2str(dataMin), num2str(dataMax)};
answer=inputdlg(prompt,name,numlines,defaultanswer);
dataMin = str2double(answer{1});
dataMax = str2double(answer{2});
%%2D or 3D?
plotTypeButton = questdlg('What kind of plot?', ...
'Question', ...
'2D', '3D', '2D');
if strcmp(plotTypeButton,'3D')
threeD = true;
else
threeD = false;
end
%%loop over the windows and conditions and plot them.
if threeD
headplot('setup', eloc, 'headSplines.spl');
end
%subplot indices
subplotInds = 1:(nWindows+2)*nCells;
subplotInds = reshape(subplotInds,nWindows+2,nCells)';
figure; hold on;
for win = 1:nWindows
for cel = 1:nCells
subplot(nCells,nWindows+2,subplotInds(cel,win+1));
if threeD
headplot(data(:,win,cel), 'headSplines.spl', ...
'maplimits', [dataMin, dataMax], ...
'electrodes', 'off', ...
'view', [azimuth, elevation]);
else
topoplot(data(:,win,cel), eloc, ...
'maplimits', [dataMin, dataMax], ...
'shading', 'interp', ...
'numcontour', 0);
end
if cel == 1
title([num2str(windows(win)) '-' num2str(windows(win+1)) 'ms']);
end
end
end
%add labels to the rows
for cel = 1:nCells
subplot(nCells,nWindows+2,subplotInds(cel,1));
set(gca,'Visible','off');
text(.5,.5,cellNames(cel),'HorizontalAlignment','Center');
end
%add colorbar
subplot(nCells,nWindows+2,sub2ind([nWindows+2,nCells],nWindows+2,1));
set(gca,'Visible','off');
h = colorbar;
set(h,'YLim',[dataMin dataMax]);
set(h,'CLim',[dataMin dataMax]);
caxis([dataMin dataMax]);
set(h,'YTick',[dataMin dataMax]);
hold off;
Memory & Cognition Blog 