Search code examples
matlabthresholdsignificance

Use findpeaks function on a time series only if this time series eclipses another


I have a question regarding findpeaks. I want to use it to detect peaks in my signal time series (Signal 1). This works fine, but I also have surrogate data, serving as a threshold of significance, of equal length (Signal 2). I now want to use findpeaks on Signal 1, but only if Signal 1 is greater then Signal 2 at that timepoint. I tried to use the regular properties of findpeaks but nothing worked so far...Here is what I have right now:

GPDC is a 9x9x512 double. Dim 1 contains partial directed coherence values estimated through a multi-variate autoregressive model in direction xi - xj, Dim 2 contains the same for xj -xi and Dim 3 represents the number of frequency bins. eEPDCsth is a 9x9x512 double containing the corresponding surrogate data. f is a 1x512 double containing the frequency values. I think right now, the >= reference isn't working, because it is not time-specific, i.e. it does not compare the signal point by point but rather in total. This is my main problem I think...

Sz=9;
for i=1:Sz
    for j=1:Sz
    if squeeze(GPDC(i,j,:)) >= squeeze(eEPDCsth(i,j,:))
       [pks_1{i,j},locs_1{i,j}] = findpeaks(squeeze(GPDC(i,j,:)),f,'npeaks',5,'MinPeakHeight', .1);
    end
    end
end

Solution

  • Here is an example that should accomplish what you have described. You didn't specify the actual contents of the 'f' vector, so I've set it to 1:512 for this example

    % data for testing
    GPDC = rand(9,9,512);
    eEPDCsth = rand(9,9,512);
    f = 1:512; % the value of the 'f' vector wasn't specified in question
    
    Sz=9;
    for i=1:Sz
        for j=1:Sz
            % find the 'raw' peaks below thresholding
            [peak_val_raw, peak_indices_raw] = findpeaks(squeeze(GPDC(i,j,:)),'npeaks',5,'MinPeakHeight', .1);
    
            % only keep peaks that are above the corresponding threshold value
            peaks_above_threshold = squeeze(GPDC(i,j,peak_indices_raw)) > squeeze(eEPDCsth(i,j,peak_indices_raw));
            peak_values_thresholded = peak_val_raw(peaks_above_threshold);
            peak_indices_thresholded = peak_indices_raw(peaks_above_threshold);
    
            pks_1{i,j} = peak_values_thresholded;
            % index into 'f' vector to match code in original question
            locs_1{i,j} = f(peak_indices_thresholded); 
    
        end
    end