Calculate duration of Peaks using findpeaks function

18 views (last 30 days)
Dear reader,
I am currently working on an analysis of storm events. I am interested in finding the duration of storms defined by windspeeds that exceed a particular threshold:
Currently, I have used the findpeaks() function to identify the location (timestamps) and maximum intensity (peak value) of the windspeeds. I have done this as follows:
[PKS,LOCS] = findpeaks(Windspeed,time,'MinPeakHeight',th)
Here, windspeed is a vector containing hourly windspeed measurements, time is vector with corresponding datenum values for the measurements of the windspeed. Finally, th is the thresholds for when windspeeds are considered a storm events.
Now, I also want to know the duration of the storm events. That is, the time between when the windspeed first exceeds th until it agains drops below th.
I have tried using 'WidthReference', however, this only gives the options 'halfprom' and 'halfheight' and both don't give the values I am looking for.
Can anyone help me with this. Many thanks in advance!

Accepted Answer

Star Strider
Star Strider on 10 Feb 2022
Try something like this —
t = linspace(0, 24*30-1, 2*30); % Simulate Time
wv = 2 + sin(2*pi*t/30) .* cos(2*pi*t/24); % Simulate Data
th = 2.4; % Set Threshold
idx = find(diff(sign(wv - th))); % Approximate Indices Where 'wv' Crosses 'th'
for k = 1:numel(idx)
idxrng = max([1,idx(k)-1]) : min([numel(t),idx(k)+1]); % Index Range For Interpolation
t_th(k,:) = interp1(wv(idxrng), t(idxrng), th); % More Precise Time Estimates For 'th' Crossings
end
% t_th
t_th = t_th(1:2*floor(numel(t_th)/2)); % Delete Incomplete Events
t_thm = reshape(t_th, 2, []); % Matrix Of Paired Events
EventDuration = diff(t_thm)
EventDuration = 1×4
45.9247 34.8986 33.6824 45.7171
figure
plot(t, wv)
hold on
plot(t, ones(size(t))*th, ':k')
plot(t_thm, ones(size(t_thm))*th, '|-r')
hold off
grid
text(mean(t_thm), ones(size(mean(t_thm)))*th, compose('Duration = %.2f\\rightarrow', EventDuration), 'Horiz','right', 'Vert','bottom', 'Rotation',-60)
It may be necessary for you to adapt this to your data. It may work as written, however it should be straightforward to change it if required.
.
  4 Comments
Michiel Smit
Michiel Smit on 11 Feb 2022
I noticed a mistake in the script above. If wv(idxrng) has succesive values that are the same (very rare but it may happen), interp1() will not work. I have solved it by adding the following lines to the for-loop for k:
% Find the threshold crossings for time and wind velocity
cross_wv = wv(idxrng);
cross_t = time(idxrng);
% Make sure all points are unique for interp1 to function
[~,id_unique] = unique(cross_wv);
% Now Approximate crossing time by running interp1()
t_th(kk,:) = interp1(cross_wv(id_unique), cross_t(id_unique), th);
Then proceed as normally as shown in the example script above.
It does not give the most accurate results, but for my work, accurate to the nearest 4 hours is more than sufficient.
Star Strider
Star Strider on 11 Feb 2022
The way I solve that is usually to interpolate to a higher resolution (more data points).
In a timetable, that can be done with the retime function, using the 'pchip' method (my favourite, unless the data are truly linear), otherwise using interp1 (again with 'pchip'). This has the disadvantage of creating data where none previously existed, however here, where the storm durations are the desired result and the time tolerance is forgiving of such inaccuracies, it likely would be appropriate and straightforward.
The problem with unique is that it elimiinates some of the data. That likely creates more problems than it solves.

Sign in to comment.

More Answers (1)

David Hill
David Hill on 10 Feb 2022
f=num2str(Windspeed>=th);
f=f(f~=' ');
[startIndex,endIndex]=regexp(f,'1*');
Duration=diff(datenum([time(startIndex),time(endIndex)])');

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!