# How to find duration of peaks/valleys in chart

8 views (last 30 days)
Sherzaad Dinah on 27 Nov 2023
Commented: Mathieu NOE on 11 Dec 2023
I have the following data from an accelerometer. I would like to find the duration of the 'main' deccelaration peaks/valleys.
I know of the 'findpeaks' function but really strugging with using it to achieve the above given my limited knowledge in matlab
would appreciate if someone could help out here.
Dyuman Joshi on 27 Nov 2023
What does 'duration' mean in this context? Please specify.
Peter Perkins on 27 Nov 2023
Your data are very noisy. You are going to need to decide what you mean by "main", and almost certainly use some kind of smoothing to find those "main" peaks/valleys.

Star Strider on 27 Nov 2023
Edited: Star Strider on 27 Nov 2023
Try this —
T1 = 180224×2 table
x[s] y,ch1[m/s^2] ______ ____________ 4 0.01202 4 -0.0206 4.0001 0.01202 4.0001 -0.01373 4.0001 0.03262 4.0001 -0.0618 4.0002 -0.09442 4.0002 -0.0206 4.0002 -0.0618 4.0003 -0.04807 4.0003 0 4.0003 0.01202 4.0004 -0.00687 4.0004 -0.00687 4.0004 -0.07382 4.0005 0.00172
VN = T1.Properties.VariableNames;
x = T1{:,1};
y = T1{:,2};
Checkx = [mean(diff(x)) std(diff(x))] % Check for Uniform Sampling
Checkx = 1×2
1.0e-04 * 0.3052 0.0222
Fs = 1/mean(diff(x))
Fs = 3.2768e+04
[yr,xr] = resample(y, x, Fs);
figure
plot(x, y)
grid
xlabel('Time (s)')
ylabel('Amplitude)')
title(["Original ‘"+string(VN{2})+"’ Signal"])
[FTyr, Fv] = FFT1(yr,xr);
figure
plot(Fv, abs(FTyr)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title(["Fourier Transform Of Resampled ‘"+string(VN{2})+"’ Signal"])
xlim([0 0.25]*1E+4)
% filt_yr = lowpass(yr, 50, Fs)
filt_yr = sgolayfilt(yr, 3, 601);
[pks,plocs] = findpeaks(filt_yr, 'MinPeakProminence',0.05);
xr_pkx = xr(plocs);
[vys,vlocs] = findpeaks(-filt_yr, 'MinPeakProminence',5);
xr_vys = xr(vlocs);
for k = 1:numel(vlocs)
vly_idx(k,2) = vlocs(k);
vly_idx(k,1) = max(plocs(plocs<vlocs(k)));
vly_idx(k,3) = min(plocs(plocs>vlocs(k)));
end
% vly_idx
figure
hp{1} = plot(x, y, 'DisplayName','Original Signal');
hold on
hp{2} = plot(x, filt_yr, 'DisplayName','Resampled & S-G Filtered Signal');
hp{3} = plot(xr(vly_idx), filt_yr(vly_idx), 'rs', 'MarkerSize',7.5, 'MarkerFaceColor','r', 'DisplayName','Valley & Border Markers');
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude)')
title(["Resampled & Filtered ‘"+string(VN{2})+"’ Signal & Analysis"])
legend([hp{1} hp{2} hp{3}(1)], 'Location','best')
xlim([5.0 7.5])
valley_xr = xr(vly_idx);
Valley_Data = array2table(valley_xr, 'VariableNames',{'Start','Minimum','End'});
Descend_xr = valley_xr(:,2) - valley_xr(:,1);
Ascend_xr = valley_xr(:,3) - valley_xr(:,2);
Total_xr = valley_xr(:,3) - valley_xr(:,1);
Valley_Data = 4×6 table
Start Minimum End Descend Ascend Total ______ _______ ______ ________ ________ ________ 5.4557 5.4802 5.504 0.024414 0.023804 0.048218 5.9696 5.9983 6.0246 0.028687 0.026276 0.054962 6.3375 6.3712 6.4029 0.03363 0.031738 0.065369 6.624 6.6514 6.6835 0.027435 0.032074 0.059509
function [FTs1,Fv] = FFT1(s,t)
s = s(:);
t = t(:);
L = numel(t);
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv);
end
The data are noisy, so the first step is to determine if the noise is broadband or band-limited. A frequency-selective filter can eliminate band-limited noise, however your data have broadband noise (as demonstrated in the Fourier transform), so either wavelet denoising or the Savitzky-Golay filter (the sgolayfilt function) is necessary to minimise it. (The sampling intervals are not constant, and actually vary significantly, so I first used the resample function to resample them to a uniform sampling frequency before calculating the Fourier transform and filtering the data.) After that, the code uses findpeaks to first locate the valleys and then the nearest peaks on either side of the valley. This is relatively straightforward, and those results are returned in the ‘vly_idx’ array. All the information in ‘Valley_Data’ are in units of ‘x’ (seconds). There are four identified peaks here. If you wantto identify more, decrease the 'MinPeakProminence' value in the second findpeaks call. It would be best to leave the rest of the code as it is, unless other data sets require that it be changed.
EDIT — Expanded the description of how the code works, and the approach I took.
.
EDIT — (27 Nov 2023 at 14:38)
My code is longer because I took time to analyse the signal first. The code itself — without the initial analysis — is simply this —
VN = T1.Properties.VariableNames;
x = T1{:,1};
y = T1{:,2};
Fs = 1/mean(diff(x))
[yr,xr] = resample(y, x, Fs);
filt_yr = sgolayfilt(yr, 3, 601);
[pks,plocs] = findpeaks(filt_yr, 'MinPeakProminence',0.05);
xr_pkx = xr(plocs);
[vys,vlocs] = findpeaks(-filt_yr, 'MinPeakProminence',5);
xr_vys = xr(vlocs);
for k = 1:numel(vlocs)
vly_idx(k,2) = vlocs(k);
vly_idx(k,1) = max(plocs(plocs<vlocs(k)));
vly_idx(k,3) = min(plocs(plocs>vlocs(k)));
end
% vly_idx
figure
hp{1} = plot(x, y, 'DisplayName','Original Signal');
hold on
hp{2} = plot(x, filt_yr, 'DisplayName','Resampled & S-G Filtered Signal');
hp{3} = plot(xr(vly_idx), filt_yr(vly_idx), 'rs', 'MarkerSize',7.5, 'MarkerFaceColor','r', 'DisplayName','Valley & Border Markers');
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude)')
title(["Resampled & Filtered ‘"+string(VN{2})+"’ Signal & Analysis"])
legend([hp{1} hp{2} hp{3}(1)], 'Location','best')
xlim([5.0 7.5])
valley_xr = xr(vly_idx);
Valley_Data = array2table(valley_xr, 'VariableNames',{'Start','Minimum','End'});
Descend_xr = valley_xr(:,2) - valley_xr(:,1);
Ascend_xr = valley_xr(:,3) - valley_xr(:,2);
Total_xr = valley_xr(:,3) - valley_xr(:,1);
If the plots are not necessary, simply remove them or comment them out so that the plots do not execute.
.

Mathieu NOE on 27 Nov 2023
hello
I have nothing against findpeaks but I prefer sometimes other solutions that are simpler and faster - therefore I usually rely on peakseek that can be found here PeakSeek - File Exchange - MATLAB Central (mathworks.com)
(I also attached here for you if it's more convenient)
main code is splitted in two sections , one with raw data and one with smoothed data (as there is indeed some noise that makes your data a bit hairy and complicates the task of find the optimal parameters for findeaks / peakseek
so my preference goes for the solution with some smoothing (bottom plot) :
clc
clearvars
t = data(:,1);
y = data(:,2);
dt = mean(diff(t));
Fs = 1/dt;
%% try with peakseek
% min time interval between peaks
time_delta_between_peaks = 0.15; % in second
minpeakdist = round(time_delta_between_peaks*Fs); % in samples
minpeakh = 0.4;
% positive peaks
[locs1, y_pks1]=peakseek(y,minpeakdist,minpeakh);
t_pks1 = t(locs1);
% negative peaks
[locs2, y_pks2]=peakseek(-y,minpeakdist,minpeakh);
t_pks2 = t(locs2);
y_pks2 = -y_pks2;
figure(1)
plot(t,y,'b',t_pks1,y_pks1,'dr',t_pks2,y_pks2,'dg');
%% do the same after some smoothing
y = smoothdata(y,'gaussian',3000);
% positive peaks
[locs1, y_pks1]=peakseek(y,minpeakdist,minpeakh);
t_pks1 = t(locs1);
% negative peaks
[locs2, y_pks2]=peakseek(-y,minpeakdist,minpeakh);
t_pks2 = t(locs2);
y_pks2 = -y_pks2;
figure(2)
plot(t,y,'b',t_pks1,y_pks1,'dr',t_pks2,y_pks2,'dg');
if you need to find the average period (in seconds) between peaks :
% average period between positive peaks
average_period_between_positive_peaks = mean(diff(t_pks1))
% average period between negative peaks
average_period_between_negative_peaks = mean(diff(t_pks2))
and you get :
average_period_between_positive_peaks = 0.3548
average_period_between_negative_peaks = 0.3448
##### 3 CommentsShow 1 older commentHide 1 older comment
Mathieu NOE on 27 Nov 2023
hello again
I am not 100% sure what you mean by duration here
this ? time delta between red & green dots
peak_duration =
0.0822
0.0781
0.0732
0.0661
0.0557
0.0445
0.0284
t = data(:,1);
y = data(:,2);
dt = mean(diff(t));
Fs = 1/dt;
% some smoothing
y = smoothdata(y,'gaussian',3000);
% find zero crossing points for negative peaks
threshold = -0.25;
[ZxP,ZxN] = find_zc(t,y,threshold)
% positive slope ZC points
t_pks1 = ZxP;
y_pks1 = interp1(t,y,t_pks1);
% negative slope ZC points
t_pks2 = ZxN;
y_pks2 = interp1(t,y,t_pks2);
% peak duration
peak_duration = t_pks1 - t_pks2;
figure(2)
plot(t,y,'b',t_pks1,y_pks1,'dr',t_pks2,y_pks2,'dg');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
Mathieu NOE on 11 Dec 2023
problem solved ?

### Categories

Find more on Multirate Signal Processing in Help Center and File Exchange

R2016b

### Community Treasure Hunt

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

Start Hunting!