Finding Peaks after Event (in a time frame)
Show older comments
I'm trying to find peaks of a dataset that occur after an event. I've attached a figure to help me explain. I am trying to find the local maximum (red arrows for examples) right after each event (marked by the dotted lines). I think I can use findpeaks for this, but I'm not sure how to set the timeframe. I want to have the following data recorded: [Event label, event time, peak height, time from event to first peak, time of peak (should be the index of the peak)], and I need the electrode label recorded.
The files I actually want to use this on are much larger (like the picture I attached), but even the zip file of one of the datasets was too big to attach, so I'm using this dataset as an example.
%% Importing_and_Plotting_EEG_Example
close all; clear all; clc;
% look to 300 ms for peak after stimulation
% stops Matlab from rounding numbers
format long;
% Importing csv data
T=readtable('Try 4.csv'); %This reads the csv data into a table named T
timestampOrig=T.Timestamp; %this makes timestampOrig equal to the original Timestamp
T.Timestamp=T.Timestamp-T.Timestamp(1); %This changes the official Timestamp column so it is 0 at the start of the recording
% Changes the sensor names from EEG_<sensor> to <sensor>
T.Properties.VariableNames(2:end)=extractAfter(T.Properties.VariableNames(2:end),'_');
%T5=head(T, 5) %shows first 5 rows of table T
Int=readtable('Event Marker Try 4_EPOCX_126368_2021.06.04T15.42.31.04.00_intervalMarker.csv');
EvtMar=Int.latency;
EvtMarLab=Int.type;
EvtLab=readtable('Try 4 Event Labels.xlsx');
table(EvtLab)
If possible, I'd like to use a for loop, similar to the one shown below, so I automate the process of finding the values for all 14 electrodes in certain sections, such as in EO, EC, CE, OM.
for i=2:width(T) %1 plot for each sensor
%finds the peaks for each sensor. The peak is currently defined as as
%point that is at least 100 uV above the values directly to the left
%and right of the peak value; this also plots the peaks
findpeaks(T{:,i},T.Timestamp, 'MinPeakProminence', 100)
xline(EvtMar,'--',{'EO','BE','BE','BE','EC','BE','CE','CE','BE','OM','OM'});
%creates ylabel as the sensor names, in bold font
ylabel(T.Properties.VariableNames(i), 'FontWeight', 'bold')
end %ends for loop
I have this section of code that I used to find the time between peaks, but it isn't quite what I'm trying to do.
%% Time Between Peaks
%Assigns the peak values and indices of peak values to variables
%<sensor>pk and <sensor>index, respectively
[AF3pk, AF3index]=findpeaks(T.AF3, T.Timestamp,'MinPeakProminence',mnpkprom); [F7pk, F7index]=findpeaks(T.F7, T.Timestamp,'MinPeakProminence',mnpkprom);
[F3pk, F3index]=findpeaks(T.F3, T.Timestamp,'MinPeakProminence',mnpkprom); [FC5pk, FC5index]=findpeaks(T.FC5, T.Timestamp,'MinPeakProminence',mnpkprom);
[T7pk, T7index]=findpeaks(T.T7, T.Timestamp,'MinPeakProminence',mnpkprom); [P7pk, P7index]=findpeaks(T.P7, T.Timestamp,'MinPeakProminence',mnpkprom);
[O1pk, O1index]=findpeaks(T.O1, T.Timestamp,'MinPeakProminence',mnpkprom); [O2pk, O2index]=findpeaks(T.O2, T.Timestamp,'MinPeakProminence',mnpkprom);
[P8pk, P8index]=findpeaks(T.P8, T.Timestamp,'MinPeakProminence',mnpkprom); [T8pk, T8index]=findpeaks(T.T8, T.Timestamp,'MinPeakProminence',mnpkprom);
[FC6pk, FC6index]=findpeaks(T.FC6, T.Timestamp,'MinPeakProminence',mnpkprom); [F4pk, F4index]=findpeaks(T.F4, T.Timestamp,'MinPeakProminence',mnpkprom);
[F8pk, F8index]=findpeaks(T.F8, T.Timestamp,'MinPeakProminence',mnpkprom); [AF4pk, AF4index]=findpeaks(T.AF4, T.Timestamp,'MinPeakProminence',mnpkprom);
%creates tables of <sensor>Indices and <sensor> Peak Values - different
%table for each sensor
AF3tab=table(AF3index, AF3pk, 'VariableNames', {'AF3 Index (s)', 'AF3 Peak (uV)'});
F7tab=table(F7index, F7pk, 'VariableNames', {'F7 Index (s)', 'F7 Peak (uV)'});
F3tab=table(F3index, F3pk, 'VariableNames', {'F3 Index (s)', 'F3 Peak (uV)'});
FC5tab=table(FC5index, FC5pk, 'VariableNames', {'FC5 Index (s)', 'FC5 Peak (uV)'});
T7tab=table(T7index, T7pk, 'VariableNames', {'T7 Index (s)', 'T7 Peak (uV)'});
P7tab=table(P7index, P7pk, 'VariableNames', {'P7 Index (s)', 'P7 Peak (uV)'});
O1tab=table(O1index, O1pk, 'VariableNames', {'O1 Index (s)', 'O1 Peak (uV)'});
O2tab=table(O2index, O2pk, 'VariableNames', {'O2 Index (s)', 'O2 Peak (uV)'});
P8tab=table(P8index, P8pk, 'VariableNames', {'P8 Index (s)', 'P8 Peak (uV)'});
T8tab=table(T8index, T8pk, 'VariableNames', {'T8 Index (s)', 'T8 Peak (uV)'});
FC6tab=table(FC6index, FC6pk, 'VariableNames', {'FC6 Index (s)', 'FC6 Peak (uV)'});
F4tab=table(F4index, F4pk, 'VariableNames', {'F4 Index (s)', 'F4 Peak (uV)'});
F8tab=table(F8index, F8pk, 'VariableNames', {'F8 Index (s)', 'F8 Peak (uV)'});
AF4tab=table(AF4index, AF4pk, 'VariableNames', {'AF4 Index (s)', 'AF4 Peak (uV)'});
pkpm="Minimum peak prominence is " + mnpkprom;
%Exports peak data to excel document - single sheet (Name: Try 4 Peaks
%and Indices)
filename = 'Try 4 Peaks and Indices.xlsx';
writematrix(pkpm, filename,'Sheet',1,'Range','A1');
writetable(AF3tab,filename,'Sheet',1,'Range','A2');
writetable(F7tab, filename, 'Sheet',1, 'Range', 'C2');
writetable(F3tab, filename, 'Sheet',1, 'Range', 'E2');
writetable(FC5tab, filename, 'Sheet',1, 'Range', 'G2');
writetable(T7tab, filename, 'Sheet',1, 'Range', 'I2');
writetable(P7tab, filename, 'Sheet',1, 'Range', 'K2');
writetable(O1tab, filename, 'Sheet',1, 'Range', 'M2');
writetable(O2tab, filename, 'Sheet',1, 'Range', 'O2');
writetable(P8tab, filename, 'Sheet',1, 'Range', 'Q2');
writetable(T8tab, filename, 'Sheet',1, 'Range', 'S2');
writetable(FC6tab, filename, 'Sheet',1, 'Range', 'U2');
writetable(F4tab, filename, 'Sheet',1, 'Range', 'W2');
writetable(F8tab, filename, 'Sheet',1, 'Range', 'Y2');
writetable(AF4tab, filename, 'Sheet',1, 'Range', 'AA2');
I also attached the full matlab code I am using in the zip file, along with the files you need to run the code.
Thank You!
EDIT/ADDED:
Below is an image of what I'm trying to get. I want the data in an excel file as a table. I will add the subject # every time, but I want to automate everything else in the table, since I have 15 datasets of about 20 mins length each to go through.
Accepted Answer
Kevin Holly
on 29 Oct 2021
Edited: Kevin Holly
on 29 Oct 2021
% Importing_and_Plotting_EEG_Example
close all; clear all; clc;
format long;
% Importing csv data
T=readtable('Try 4.csv'); %This reads the csv data into a table named T
timestampOrig=T.Timestamp; %this makes timestampOrig equal to the original Timestamp
T.Timestamp=T.Timestamp-T.Timestamp(1); %This changes the official Timestamp column so it is 0 at the start of the recording
% Changes the sensor names from EEG_<sensor> to <sensor>
T.Properties.VariableNames(2:end)=extractAfter(T.Properties.VariableNames(2:end),'_');
%T5=head(T, 5) %shows first 5 rows of table T
Int=readtable('Event Marker Try 4_EPOCX_126368_2021.06.04T15.42.31.04.00_intervalMarker.csv');
EventTime=Int.latency;
EventMarkerLabel=Int.type;
EvtLab=readtable('Try 4 Event Labels.xlsx');
EvtLab{2,2} = {'Blink Eye'};
table(EvtLab)
% Preallocate variable
Event_Labels = cell(length(EventTime),1);
% Identify Event Labels
for i = 1:size(EvtLab,1)
Event_Labels(contains(EventMarkerLabel(1:length(EventTime)),EvtLab{i,2},'IgnoreCase',true),1) = EvtLab{i,1};
end
Plots all together
p=plot(T.Timestamp,T{:,2:width(T)});
xline(EventTime,'--',Event_Labels);
% Preallocate variables
Firstpeakafterevent = zeros([width(T), size(EventTime)]);
timeoffirstpeakafterevent = zeros([width(T), size(EventTime)]);
timefromeventtofirstpeak = zeros([width(T), size(EventTime)]);
% Plot
for i=2:width(T)
[pks,locs] = findpeaks(T{:,i},T.Timestamp, 'MinPeakProminence', 100);
for ii = 1:length(EventTime)
index = find(locs>EventTime(ii),1);
if ~isempty(index)
Firstpeakafterevent(i-1,ii,:) = pks(index);
timeoffirstpeakafterevent(i-1,ii,:) = locs(index);
timefromeventtofirstpeak(i-1,ii,:) = locs(index)-EventTime(ii);%time from event to first peak
end
end
x = timeoffirstpeakafterevent(i-1,:);
y = Firstpeakafterevent(i-1,:);
x(x==0) = [];
y(y==0) = [];
hold on
scatter(x,y,'v','filled','MarkerFaceColor',p(i-1).Color)
end
legend(p,T.Properties.VariableNames(2:width(T)))
If you want separate figures
for i=2:width(T)
[pks,locs] = findpeaks(T{:,i},T.Timestamp, 'MinPeakProminence', 100);
for ii = 1:length(EventTime)
index = find(locs>EventTime(ii),1);
if ~isempty(index)
Firstpeakafterevent(i-1,ii,:) = pks(index);
timeoffirstpeakafterevent(i-1,ii,:) = locs(index);
timefromeventtofirstpeak(i-1,ii,:) = locs(index)-EventTime(ii);%time from event to first peak
end
end
x = timeoffirstpeakafterevent(i-1,:);
y = Firstpeakafterevent(i-1,:);
x(x==0) = [];
y(y==0) = [];
figure(i-1)
plot(T.Timestamp,T{:,i})
hold on
scatter(x,y,'v','filled')
xline(EventTime,'--',Event_Labels);
ylabel(T.Properties.VariableNames(i))
end
11 Comments
Christina Diersing
on 2 Nov 2021
Actually, I was trying to get an excel file/table like the one below. I can add the subject number, but I'm trying to get the rest of the table with an automated process, since I have to go through 15 datasets of about 20 minutes each.
Christina Diersing
on 3 Nov 2021
Kevin Holly
on 5 Nov 2021
Edited: Kevin Holly
on 5 Nov 2021
You can create a table of the data as such:
timefromeventtofirstpeak = timefromeventtofirstpeak'; %transposing data so it fits table
timeoffirstpeakafterevent = timeoffirstpeakafterevent';
Firstpeakafterevent = Firstpeakafterevent';
t=table
t.EventLabels = Event_Labels
t.Latency = timefromeventtofirstpeak(:,1);% you can change 1 into a variable in used for a for loop
t.Time_of_Peak = timeoffirstpeakafterevent(:,1);
t.PeakHeight = Firstpeakafterevent;
You can save this in Excel using the writetable function. You could create a for loop to create a table for each electrode channel. You could either combine the tables or save each table on a separate Excel sheet.
Christina Diersing
on 9 Nov 2021
Kevin Holly
on 9 Nov 2021
I preallocated the variables incorrectly.
It should be "width(T)-1" instead of "width(T)" since I started the for loop at 2 instead of 1 because I believe the first column was header information.
Firstpeakafterevent = zeros([width(T)-1, size(EventTime)]);
timeoffirstpeakafterevent = zeros([width(T)-1, size(EventTime)]);
timefromeventtofirstpeak = zeros([width(T)-1, size(EventTime)]);
Christina Diersing
on 9 Nov 2021
Kevin Holly
on 9 Nov 2021
Edited: Kevin Holly
on 9 Nov 2021
check the size of Firstpeakafterevent
size(Firstpeakafterevent)
Does the amount of variable names match the columns? The output of the command aboves gives the numbers of rows then the number of columns.
You may need to place that command before it was transposed.
pkHeighttab=table(Firstpeakafterevent, 'VariableNames',{'AF3 Pk','F7 Pk','F3 Pk','FC5 Pk','T7 Pk','P7 Pk','O1 Pk','O2 Pk','P8 Pk','T8 Pk','FC6 Pk','F4 Pk','F8 Pk','AF4 Pk'});
Firstpeakafterevent = Firstpeakafterevent';
Why are you trying to create this new table?
Christina Diersing
on 9 Nov 2021
The way you created the table, it placed the double array under one column.
Firstpeakafterevent=rand(11,14);
pkHeighttab=table(Firstpeakafterevent)
You would need to separate each column by a comma as such:
pkHeighttab=table(Firstpeakafterevent(:,1),Firstpeakafterevent(:,2),Firstpeakafterevent(:,3))
Or you could convert the double array to a table
pkHeighttab = array2table(Firstpeakafterevent)
pkHeighttab.Properties.VariableNames = {'AF3 Pk','F7 Pk','F3 Pk','FC5 Pk','T7 Pk','P7 Pk','O1 Pk','O2 Pk','P8 Pk','T8 Pk','FC6 Pk','F4 Pk','F8 Pk','AF4 Pk'}
Kevin Holly
on 9 Nov 2021
For the electrodes, I was envisioning that you would use a for loop to create a table for each electrode, where each table would be saved to its own sheet. Where the electrode label would be saved in cell 'B4' and on the sheet name itself for instance.
Something like this:
filename = 'test.xlsx';
Electrodes = {'AF3','F7','F3','FC5','T7','P7','O1','O2','P8','T8','FC6','F4','F8','AF4'};
for i=1:14
%insert code to generate table%
writetable(t,filename,'Sheet',Electrodes{i},'Range','A5');
writematrix("Electrodes",filename,'Sheet',Electrodes{i},'Range','A4')
writematrix(Electrodes{i},filename,'Sheet',Electrodes{i},'Range','B4')
writematrix("Subject #",filename,'Sheet',Electrodes{i},'Range','A4')
end
Christina Diersing
on 10 Nov 2021
More Answers (0)
Categories
Find more on EEG/MEG/ECoG in Help Center and File Exchange
on 29 Oct 2021
on 10 Nov 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
