Time to peak using findpeaks

I'm using findpeaks to locate multiple peaks in my function but I want to know the rising time to the peak.

Findpeaks gives you the 'width' output but its not working because it suposses that the peak is in the middle of the valleys (and its based on the prominece, it would be nice to be referenced to lowestPoint/left valley) but most of these peaks are not symetric so its not working properly, any ideas?

Answer 1

[Star Strider]

I usually have findpeaks return the indices of the peaks.

Your definition of the start time point isnot obvious. I would use interp1 to define it as a function of its amplitude, then subtract the time returned from the time of the peak.

Example --

t = datetime('07:00:00', InputFormat='HH:mm:ss') : minutes(1) : datetime('08:00:00', InputFormat='HH:mm:ss');
t.Format = 'HH:mm:ss';
t = t(:);

mp = round(numel(t)/2);

f = exp(-minutes(t - t(mp)).^2/50) * 1000 + 2500;           % Create Data (Curve)

[pks,locs] = findpeaks(f)                                   % Peak & Index

pks = 3500

locs = 31

thrshld = 2550;                                             % Amplitude Threshold
idxrng = 1 : locs;                                          % Index Rnage For Interpolation
thrshld_time = interp1(f(idxrng), t(idxrng), thrshld)       % Interpolate To Find Associated Time

thrshld_time = datetime
   07:17:43

t_interval = t(locs) - thrshld_time

t_interval = duration
   00:12:16

figure
plot(t, f)
hold on
plot([thrshld_time thrshld_time], [2550 2750], ":r")
plot([thrshld_time t(locs)], 2750*[1 1], '|-r')
hold off
grid
yline(thrshld, '--r')
xline(t(locs), '--r')
text(thrshld_time+t_interval/2, 2750, sprintf('\\leftarrow \\Delta = %s', t_interval), Rotation=45)
axis('padded')

●

.

[David Santos]

thanks for your answer!

It is a good idea the problem is that the threshold value is no "predictable", I mean, is something that changes a lot between each peak of my function

[Star Strider]

My pleasure!

It would help to have your data, as well as a description of the problem you want to solve, and whatever characteristics define the threshold value. I can probably include code for these, if I understand how you want to define them.

[David Santos]

Great! T

I've attached a table with some data example.

I'm interested in the peaks information, which relates to the closing of a mussel valve ("mouth"). I want to know the speed of each of these closings which can be seen as the time from the "base level" to the peak.

This 'base level' changes over the time, even more than the data attached but I would be happy with a good approximation.. I think that the characteristic that defines it is the point before the peak where the slope is close to zero. In the next figure I've marked some of the base levels so you can figure want I mean for the rest of the peaks.

Zooming (3rd peak):

I've tried with risetime function but if fails, i think because sometimes the rise is not a straight line but I'm not sure... but the way these function approaches the problem seems good to me. Maybe something in this direction?

All the best

[Star Strider]

Thank you!

I don't see the table (I assume an Excel or .csv file). It would be best to have the time vector as well as the signal vector.

The first part of analysing your data would be to remove the baseline drift and baseline offset. That would require using the Signal Processing Toolbox highpass function. (If you do not have that function, I can easily write equivalent code that should work in earlier releases.)

After that, determining the rise time should be relatively straightforward. Once the start time and the peak time are determined for each pulse, those can be translated back to the original data if necessary, and those amplitude values determined. It would also be straightforward to return the values for the filtered data.

You could also use findpeaks and the Statistics and Machine Learning Toolbox functions to get the Poisson distribution parameters of the peaks, if desired.

I will wiat for the file tto appear.

.

[David Santos]

SORRY!

I forgot the attachment..

[Star Strider]

Thank you!

The filtering approach unfortunately will not work. Although it removes the baseline offset and much of the baseline drift, it distorts the signal too much.

The best alternative then is to detect the peaks and then use the gradient (numerical derivative) function to find the first significant increase in the gradient in the 90 index positions preceding the peak. (I chose the '90' value empirically. It may be necessary to adjust that for other data sets.)

That seems to work here. I did not need to do any interpolation, since part of the code involved detecting the rise in the gradient result to define the start point for each peak. The time interval then was simply the difference between the start point and the peak. (I plotted them out individually to be certain the results were reasonable. I left that in the code, commented-out.) The '-2' in the idxv calculation is a sort of 'fudge factor' to correct for the gradient result not being exact.

My revised code --

LD = load('dataExample.mat')

LD = struct with fields:
    Tv: [37490×1 timetable]

Tv = LD.Tv

Tv = 37490×1 timetable
           timedt             C00  
    _____________________    ______

    06/25/25 16:54:59.000    1754.4
    06/25/25 16:55:00.000    1840.8
    06/25/25 16:55:01.000    1833.1
    06/25/25 16:55:02.000    1835.9
    06/25/25 16:55:03.000    1838.3
    06/25/25 16:55:04.000    1837.7
    06/25/25 16:55:05.000    1835.9
    06/25/25 16:55:06.000    1835.9
    06/25/25 16:55:07.000    1831.8
    06/25/25 16:55:08.000    1835.9
    06/25/25 16:55:09.000    1837.5
    06/25/25 16:55:10.000    1831.1
    06/25/25 16:55:11.000    1838.3
    06/25/25 16:55:12.000    1831.7
    06/25/25 16:55:13.000    1835.9
    06/25/25 16:55:14.000    1832.8

timedt = Tv.timedt;
C00 = Tv.C00;

figure
plot(timedt, C00)
grid
xlabel('Time')
ylabel('Amplitude')
title('Original Signal')

figure
plot(timedt, C00)
grid
xlabel('Time')
ylabel('Amplitude')
title('Original Signal')
xlim([min(timedt) datetime(2025,06,25,20,00,00)])

t_sec = seconds(timedt - timedt(1));
Fs = 1/mean(diff(t_sec))

Fs = 1

Fsd = std(1/mean(diff(t_sec)))

Fsd = 0

[FTs1,Fv] = FFT1(C00,t_sec);

Fco = 3E-4;

figure
semilogx(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform')
xlim([0 0.02])
xline(Fco, '--r', 'F_{co}')

●

C00_filt = highpass(C00, Fco, Fs, ImpulseResponse='iir');

figure
plot(timedt, C00_filt)
grid
xlabel('Time')
ylabel('Amplitude')
title('Filtered Signal')

% figure
% plot(timedt, C00_filt)
% grid
% xlabel('Time')
% ylabel('Amplitude')
% title('Filtered Signal')
% xlim([min(timedt) datetime(2025,06,25,20,00,00)])

[pks, locs] = findpeaks(C00, MinPeakProminence=150);
NrPks = numel(locs)

NrPks = 22

% figure
% tiledlayout(6,4)
for k = 1:numel(locs)
    % k    
    idxrng = max(1,locs(k)-90) : locs(k)+5;
    dsdt = gradient(C00(idxrng));
    idx = find(dsdt >= 10, 1, 'first');
    minv = C00(idx);
    idxv(k,:) = min(idxrng)+idx-2;
    
    t_interval(k,:) = timedt(locs(k)) - timedt(idxv(k));
    
    % figure
    % % nexttile
    % yyaxis left
    % plot(timedt(idxrng), C00(idxrng))
    % yyaxis right
    % plot(timedt(idxrng), dsdt)
    % grid
    % xline(timedt(idxv(k)),'--r')
    % title("Peak "+k)
end

figure
plot(timedt, C00, DisplayName="Data")
hold on
plot(timedt(locs), C00(locs), 'vr', MarkerFaceColor='r', DisplayName="Start")
plot(timedt(idxv), C00(idxv), '^g', MarkerFaceColor='g', DisplayName='Peak')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title('Original Signal')
legend(Location='best')

Results = table((1:numel(locs)).', timedt(idxv), timedt(locs), t_interval, seconds(t_interval), VariableNames=["Peak Nr","Start Time","Peak Time","Time Interval","Time Interval (s)"])

Results = 22×5 table
    Peak Nr         Start Time                Peak Time          Time Interval    Time Interval (s)
    _______    _____________________    _____________________    _____________    _________________

       1       06/25/25 17:56:21.000    06/25/25 17:56:44.000      00:00:23              23        
       2       06/25/25 17:56:21.000    06/25/25 17:56:51.000      00:00:30              30        
       3       06/25/25 18:33:23.000    06/25/25 18:33:49.000      00:00:26              26        
       4       06/25/25 18:33:23.000    06/25/25 18:33:56.000      00:00:33              33        
       5       06/25/25 19:09:47.000    06/25/25 19:11:06.000      00:01:19              79        
       6       06/25/25 20:27:08.000    06/25/25 20:27:31.000      00:00:23              23        
       7       06/25/25 20:44:06.000    06/25/25 20:44:26.000      00:00:20              20        
       8       06/25/25 21:57:34.000    06/25/25 21:59:05.000      00:01:31              91        
       9       06/25/25 22:45:31.000    06/25/25 22:46:04.000      00:00:33              33        
      10       06/25/25 23:03:12.000    06/25/25 23:03:42.000      00:00:30              30        
      11       06/25/25 23:26:53.000    06/25/25 23:27:25.000      00:00:32              32        
      12       06/25/25 23:34:35.000    06/25/25 23:35:33.000      00:00:58              58        
      13       06/25/25 23:58:42.000    06/25/25 23:59:35.000      00:00:53              53        
      14       06/26/25 00:22:39.000    06/26/25 00:23:30.000      00:00:51              51        
      15       06/26/25 00:22:39.000    06/26/25 00:23:48.000      00:01:09              69        
      16       06/26/25 00:43:53.000    06/26/25 00:44:36.000      00:00:43              43        

[lambdahat,lambdahatci] = poissfit(Results.("Time Interval (s)"))

lambdahat = 39.6818

lambdahatci = 2×1

37.0495 42.3141

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function [FTs1,Fv] = FFT1(s,t)
% One-Sided Numerical Fourier Transform
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

.

[David Santos]

Thanks a lot this is a great work!

[Star Strider]

As always, my pleasure!

Thank you!