findpeaks () doesn't find all peaks in function

Hello all,

i have recently written a function that is supposed to detect when an upwards slope in my data is deviated using the derivative of that slope. A matrix of concatenated EEG data epochs is supposed to be analyzed like this and it works for the first few epochs. However, I know that in the 6st epoch there is supposed to be a deviation, but the program doesnt find this and afterwards sais that the derivative of an upwards slope doesnt have a peak, which is mathematically just not logical from my understanding. Therefore, I am suspecting that there is a shift in the data at some point, but I have spent the whole day looking and couldnt find it.

Id really appreciate any advice or input. Thank you in advance.

function [labelLocs,epochLocs,stdpeak] = ExtrN2b(x,freq,varargin)

%  x is a matrix where each column contains data corresponding to a
%  channel.

% Only applicable if epoch size is 1sec
for l = 1:size(x,2)
    num = 1;
    supercount = 1;
    megacount = 1;
    for i = 1:(size(x(:,l),1)/freq)
        [yMax,xMax] = findpeaks(x(num:(num+(freq-1)),l));
        [yMin,xMin] = findpeaks(-(x(num:(num+(freq-1)),l)));
        if xMax(1,1) < xMin(1,1)
            xMax(1) = [];
        end
        for h = 1:(size(xMax,1))
            d = diff(x(xMin(h):(xMax(h)+1),l));
            if size(d,1) >= 3
                [ydiff,xdiff] = findpeaks(d); %doesnt detect all peaks, e.g. VP01f8 epoch6, and sometimes detects 0 peaks (shift in data?)
                stdpeak(supercount,l) = size(xdiff,1);
                if size(xdiff,1) > 1
                    labelLocs(supercount,l) = 1;
                    epochLocs(megacount,l) = 1;
                else
                    labelLocs(supercount,l) = 0;
                end
            else
            end
            supercount = supercount+1;
            d(1:end) = [];
            xdiff(1:end) = [];
        end
        num = num+freq;
        megacount = megacount +1;
        xMax(1:end) = [];
        xMin(1:end) = [];
    end
end
k = size(labelLocs,2);
for y = 1:size(labelLocs,1)
    if labelLocs(y,1:k) == ones(1,k)
        labelLocs(y,k+1) = 1;
    else
        labelLocs(y,k+1) = 0;
    end
end
end

[Dyuman Joshi]

Which inputs are used for calling the function?

What is the relevance/context of the attached .mat file?

[Michel]

Hi, the provided data is an excerpt of the input I use

Answer 1

[William Rose]

@Michel,

I see that the mat file contains the vector VP01f8, size 28600x1. I assume this is the signal of interest. I will refer to it as x(t). Since you have not specified a sampling rate, I assume fs=1 kHz.

load VP01f8; x=VP01f8; fs=1000; t=(0:length(x)-1)/fs;
plot(t,x,'-b'); grid on; xlabel('Time (s)'); ylabel('x(t)')

x(t) has many sharp peaks and valleys. findpeaks(x), with default settings, finds 975 peaks. findpeaks(-x) finds 974 peaks. What exactly are you looking for, when you have such a large number of peaks?

Explain the outputs from your function. What is the meaning of the figure it generates (which has no labels).

[labelLocs,epochLocs,stdpeak]=ExtrN2b(x,fs);

figure; subplot(311); plot(labelLocs(:,1),'-r.'); title('labelLocs(1,:)');
subplot(312); plot(labelLocs(:,2),'-g.'); title('labelLocs(2,:)');
subplot(313); plot(stdpeak,'-b.'); title('stdpeak');

What do these outputs mean? What do they tell you about the EEG?

The un-plotted output, epochLocs, is a vector of 18 zeros and a one. What is it for?

You write: "...a function to detect when an upwards slope in my data is deviated using the derivative of that slope" Is your real interrest in peaks in dx/dt? Or in d2x/dt2 (the derivative of the slope)? Your code uses findpeaks() to find peaks in x(t) and in -x(t), and, later, in dx/dt. Your code does not use findpeaks() to detect peaks in d2x/dt2. Please clarify what you mean by "using the derivative of that slope".

Please explain why you want to find the peaks, since that may lead to a useful discussion of how to accomplish your ultimate goal.

You write: "the program ... sais that the derivative of an upwards slope doesnt have a peak, which is mathematically just not logical from my understanding". The code does not say this, directly or indirectly. It never looks for peaks in d2x/dt2.

You write "I am suspecting that there is a shift in the data at some point, but I have spent the whole day looking and couldnt find it". I do not understand why you suspect there is a shift in the data. The plot of x(t), above, shows no shift.

Style suggestions:

1. Lower case L, which you use as a loop variable and channel number index, looks like the numeral 1, in the default code font, which makes it easier to make mistakes (ask me how I know), and it makes it harder to understand and debug code. Use something else.
2. The code includes

[yMax,xMax] = findpeaks(x(num:(num+(freq-1)),l));

and similar lines. This is confusing, because ymax contains the values of the peaks of x, and xmax contains the locations of the peaks in x. I think

[xMax,tMax] = findpeaks(x(num:(num+(freq-1)),l),fs);

or

[xMax,kMax] = findpeaks(x(num:(num+(freq-1)),l));

would make your code easier to understand.

3. Delete vargin, since it is never used.

[Michel]

Thank you very much for your help first of all. I really appreciate it. I apologize for the confusion, I wanted to keep my description brief, but I‘ll do my best to explain my intentions.

The purpose of this program is to detect an event related potential (ERP) which has a slight „bump“ in the upwards slope as it’s defining feature. It‘s a N2b/P3 complex if that means something to you.

To achieve this the code first divides the data back into the single epochs (1sec long data measurements after an impulse) and computes the positive and negative peaks. Then the upwards slopes are looked at. For each upwards slope the first derivation is calculated, assuming that if there is a deviation in the upwards slope the first derivation would exhibit two peaks. Examples for this is can be found in the first, third and sixth epoch.

labelLocs then notes is an upwards slope meets these criteria Stdpeaks notes the number of peaks detected in the first derivation

Stdpeaks will be used as a predictor to train a linear discriminate analysis algorithm in the next step labelLocs will be used as labels to train the algorithm

The input data I provided is a participants EEG recording from channel f8. Sampling frequency was 200Hz.

I hope this cleared some things up. And again, I am very grateful for your help. If there are any other questions, please ask and I’ll do my best to explain my intentions.

[William Rose]

@Michel,

Since fs=200 Hz, the recording duration is 143 s. What is the timing of the triggering events? Is there just one event at t=0, and none thereafter? Or one event per second, at t=0, 1, 2, ..., for the entire recording? Or something else?

Is it your aim to estimate whether or not there was an N2b/P3 response after each stimulus? (I.e. the N2b response to each stimulus is 0 or 1.) Or do you assume there's an N2b/P3 after every stimulus, and you want to determine its amplitude?

You said the ERP of interest "has a slight „bump“ in the upwards slope as it’s defining feature". Please cite the paper that describes this defining feature, so that I can understand if your code is effectively searching for the defining feature.

N2b latency is considered to be 200-400 ms by some (e.g. Czligler et al 1996), while others consider the latency to be narrower, 245-290 ms (e.g. Senkowski & Herrmann 2002). The P3 (or P300) latency to peak is similar. See review by Polich, 2007, here. How do you want to define the N2b? What paper do you consider to be the best, or the standard in the field, for describing how to detect or quantify the N2b/P3 response?

Let's take this offline. Please contact me by email by clicking on the envelope icon that appears at the top right of the window that pops up when you click the "WR" circle next to my answers.