# How to find duration of peaks/valleys in chart

6 views (last 30 days)

Show older comments

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.

##### 2 Comments

Peter Perkins
on 27 Nov 2023

### Answers (2)

Star Strider
on 27 Nov 2023

Edited: Star Strider
on 27 Nov 2023

Try this —

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

VN = T1.Properties.VariableNames;

x = T1{:,1};

y = T1{:,2};

Checkx = [mean(diff(x)) std(diff(x))] % Check for Uniform Sampling

Fs = 1/mean(diff(x))

[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]

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.

.

##### 0 Comments

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

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

##### 3 Comments

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

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!