How to find duration of peaks/valleys in chart

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]

What does 'duration' mean in this context? Please specify.

[Peter Perkins]

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.

Answer 1

[Star Strider]

Try this —

T1 = readtable('NV1-9999.csv', 'VariableNamingRule','preserve')

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);
Add_Data = array2table([Descend_xr Ascend_xr Total_xr],'VariableNames',{'Descend','Ascend','Total'});

Valley_Data = [Valley_Data Add_Data]

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 —

T1 = readtable('NV1-9999.csv', 'VariableNamingRule','preserve')

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);

Add_Data = array2table([Descend_xr Ascend_xr Total_xr],'VariableNames',{'Descend','Ascend','Total'});

Valley_Data = [Valley_Data Add_Data]

If the plots are not necessary, simply remove them or comment them out so that the plots do not execute.

.

Answer 2

[Mathieu NOE]

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

data = readmatrix('NV1-9999.csv');
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

[Sherzaad Dinah]

while this solution is much simple that the previous one, if does not compute the duration for the peak.

but I guess it woud be a matter of find the zero-crossing points for those 'peaks' and caculated the difference therefrom ;)

[Mathieu NOE]

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

data = readmatrix('NV1-9999.csv');
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]

hello again @Sherzaad Dinah

problem solved ?