Time to peak using findpeaks
Show older comments
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?
Accepted Answer
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
thrshld = 2550; % Amplitude Threshold
idxrng = 1 : locs; % Index Rnage For Interpolation
thrshld_time = interp1(f(idxrng), t(idxrng), thrshld) % Interpolate To Find Associated Time
t_interval = t(locs) - thrshld_time
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')
.
8 Comments
David Santos
on 15 Sep 2025
Star Strider
on 15 Sep 2025
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
on 15 Sep 2025
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
on 15 Sep 2025
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
on 16 Sep 2025
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')
Tv = LD.Tv
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))
Fsd = std(1/mean(diff(t_sec)))
[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)
% 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)"])
[lambdahat,lambdahatci] = poissfit(Results.("Time Interval (s)"))
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
on 16 Sep 2025
Star Strider
on 16 Sep 2025
As always, my pleasure!
Thank you!
More Answers (0)
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
