eeg_peaks - Find peaks in EEG data
Usage: [p] = eeg_peaks(p)
The p structure is explained in eeg_toolbox_defaults.
For this function, it must contain the field called
'volt.data,' which is assumed to be a matrix
with N rows of EEG sample values from M columns
of electrodes.
The indices to peak values are returned to
p.volt.peaks.data (which is size(p.volt.data)).
The values are:
0 no peak
1 a +ve peak
-1 a -ve peak
This function finds all peaks in a waveform, where
f(x) changes from an increasing to a decreasing
function and vice versa). These are local minima
and maxima that have a first derivate of zero and
satisfy the first derivative test.
Thus, the +ve/-ve peak refers to peaks and troughs
in a waveform, regardless of the +ve/-ve value of
the waveform at that point. The peak values are
commonly +ve/-ve potentials, but they may not be.
A +ve peak can have a -ve value where the change
from increasing f(x) to decreasing f(x) occurs
entirely in the range of negative values.
All the +ve/-ve peak values are returned in:
p.volt.peaks.all
p.volt.peaks.pos
p.volt.peaks.neg
If the p.volt.timeArray is defined, then the timing
of the peaks is returned in:
p.volt.peaks.alltimes
p.volt.peaks.postimes
p.volt.peaks.negtimes
eg, p.volt.timeArray = linspace(-20,20,200)';
p.volt.data = sin(linspace(-20,20,200))';
p.volt.data(:,2) = sin(-1 .* linspace(-20,20,200))';
[p] = eeg_peaks(p);
scatter(p.volt.timeArray(p.volt.peaks.all),p.volt.data(p.volt.peaks.all))
hold on
plot(p.volt.timeArray,p.volt.data)0001 function [p] = eeg_peaks(p) 0002 0003 % eeg_peaks - Find peaks in EEG data 0004 % 0005 % Usage: [p] = eeg_peaks(p) 0006 % 0007 % The p structure is explained in eeg_toolbox_defaults. 0008 % For this function, it must contain the field called 0009 % 'volt.data,' which is assumed to be a matrix 0010 % with N rows of EEG sample values from M columns 0011 % of electrodes. 0012 % 0013 % The indices to peak values are returned to 0014 % p.volt.peaks.data (which is size(p.volt.data)). 0015 % The values are: 0016 % 0 no peak 0017 % 1 a +ve peak 0018 % -1 a -ve peak 0019 % 0020 % This function finds all peaks in a waveform, where 0021 % f(x) changes from an increasing to a decreasing 0022 % function and vice versa). These are local minima 0023 % and maxima that have a first derivate of zero and 0024 % satisfy the first derivative test. 0025 % 0026 % Thus, the +ve/-ve peak refers to peaks and troughs 0027 % in a waveform, regardless of the +ve/-ve value of 0028 % the waveform at that point. The peak values are 0029 % commonly +ve/-ve potentials, but they may not be. 0030 % A +ve peak can have a -ve value where the change 0031 % from increasing f(x) to decreasing f(x) occurs 0032 % entirely in the range of negative values. 0033 % 0034 % All the +ve/-ve peak values are returned in: 0035 % 0036 % p.volt.peaks.all 0037 % p.volt.peaks.pos 0038 % p.volt.peaks.neg 0039 % 0040 % If the p.volt.timeArray is defined, then the timing 0041 % of the peaks is returned in: 0042 % 0043 % p.volt.peaks.alltimes 0044 % p.volt.peaks.postimes 0045 % p.volt.peaks.negtimes 0046 % 0047 % eg, p.volt.timeArray = linspace(-20,20,200)'; 0048 % p.volt.data = sin(linspace(-20,20,200))'; 0049 % p.volt.data(:,2) = sin(-1 .* linspace(-20,20,200))'; 0050 % [p] = eeg_peaks(p); 0051 % scatter(p.volt.timeArray(p.volt.peaks.all),p.volt.data(p.volt.peaks.all)) 0052 % hold on 0053 % plot(p.volt.timeArray,p.volt.data) 0054 % 0055 0056 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:33 $ 0057 0058 % Licence: www.gnu.org GPL, no implied or express warranties 0059 % Created: 02/2002, Darren.Weber_at_radiology.ucsf.edu 0060 % - adapted a perl foreach loop to peaks_matrix 0061 % to take advantage of matlab speed. 0062 % 05/2002, Darren.Weber_at_radiology.ucsf.edu 0063 % - used matlab 'diff' command 0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0065 0066 fprintf('EEG_PEAKS\n...calculating peaks...'); tic; 0067 0068 p.volt.peaks.data = peaks_matrix(p.volt.data); 0069 0070 % Find all the peak values 0071 p.volt.peaks.all = p.volt.data .* (p.volt.peaks.data ~= 0); 0072 p.volt.peaks.pos = p.volt.data .* (p.volt.peaks.data > 0); 0073 p.volt.peaks.neg = p.volt.data .* (p.volt.peaks.data < 0); 0074 0075 if ~isfield(p.volt,'timeArray'), 0076 fprintf('\n...p doesn''t contain ''volt.timeArray''.\n'); 0077 return 0078 elseif isempty(p.volt.timeArray), 0079 fprintf('\n...p.volt.timeArray is empty.\n'); 0080 return 0081 end 0082 0083 % make sure that volt.timeArray is same size as volt.data 0084 if ~isequal(size(p.volt.timeArray),size(p.volt.data)), 0085 p.volt.timeArray = repmat(p.volt.timeArray(:,1),1,size(p.volt.data,2)); 0086 end 0087 % Find the peak timing 0088 p.volt.peaks.alltimes = p.volt.timeArray .* (p.volt.peaks.data ~= 0); 0089 p.volt.peaks.postimes = p.volt.timeArray .* (p.volt.peaks.data > 0); 0090 p.volt.peaks.negtimes = p.volt.timeArray .* (p.volt.peaks.data < 0); 0091 0092 t = toc; 0093 fprintf('done (%5.2f sec)\n',t); 0094 return 0095 0096 0097 0098 function [peaks] = peaks_matrix(data) 0099 0100 % This function should be quicker than that below. 0101 % The code below remains as it has been proved to 0102 % be reliable, but this code effectively replaces it. 0103 % Both code sets produce the same results. 0104 0105 [r,c] = size(data); 0106 % replicate data, but staggered one point 0107 data2 = zeros(1,c); 0108 data2(2:r+1,:) = data; 0109 % add zero values to end of data points 0110 data(r+1,:) = zeros(1,c); 0111 % generate boolean difference between two matrices 0112 dif = data < data2; 0113 0114 % Now replicate this boolean matrix 0115 [br,bc] = size(dif); 0116 dif2 = zeros(1,bc); 0117 dif2(2:br+1,:) = dif; 0118 % add zeros to end of boolean dif 0119 dif(br+1,:) = 0; 0120 0121 % Calculate the peaks 0122 peaks = dif(2:r+1,:) - dif2(2:r+1,:); 0123 0124 % now trim peaks dif back to size of data 0125 peaks(1,:) = 0; % first point cannot be a peak 0126 peaks(r,:) = 0; % last point cannot be a peak 0127 0128 0129 0130 % Could use matlab diff command, like this 0131 %[r,c] = size(data); 0132 %decdat = (diff(data)<=0) .* -1; 0133 %incdat = diff(data)>=0; 0134 0135 %zrow = zeros(1,c); 0136 %dirdat1 = [ zrow; decdat + incdat ]; 0137 %dirdat2 = [ decdat + incdat; zrow ]; 0138 0139 %negpeaks = ((dirdat1 - dirdat2) < 0) .* -1; 0140 %pospeaks = (dirdat1 - dirdat2) > 0; 0141 %peaks = negpeaks + pospeaks; 0142 0143 %peaks(1,:) = 0; % first point cannot be a peak 0144 %peaks(end,:) = 0; % last point cannot be a peak 0145 0146 0147 return 0148 0149 0150 % The function 'peaks_matrix' above replaced the 0151 % following for loop (02/2002, Darren Weber) 0152 0153 % foreach electrode 0154 %for e = 1:size(p.volt.data,2), 0155 % v = p.volt.data(:,e); 0156 % n = 1; 0157 % volt = v(n); 0158 % while (n <= size(v,1)), 0159 % if(volt < v(n)), 0160 % volt = v(n); n = n + 1; 0161 % if(volt < v(n)), p.volt.peaks.data(n-1,e) = 0; 0162 % else p.volt.peaks.data(n-1,e) = 1; 0163 % end 0164 % elseif(volt > v(n)), 0165 % volt = v(n); n = n + 1; 0166 % if(volt > v(n)), p.volt.peaks.data(n-1,e) = 0; 0167 % else p.volt.peaks.data(n-1,e) = -1; 0168 % end 0169 % else 0170 % volt = v(n); n = n + 1; 0171 % p.volt.peaks.data(n-1,e) = 0; 0172 % end 0173 % if isequal(n,size(v,1)), 0174 % break; 0175 % end 0176 % end 0177 % p.volt.peaks.data(1,e) = 0; % first point cannot be a peak 0178 % p.volt.peaks.data(n,e) = 0; % last point cannot be a peak 0179 %end 0180 %p.volt.peaks.all = find(p.volt.peaks.data ~= 0); 0181 %p.volt.peaks.pos = find(p.volt.peaks.data > 0); 0182 %p.volt.peaks.neg = find(p.volt.peaks.data < 0); 0183 %return