How to find duration of peaks/valleys in chart

Question

Sherzaad Dinah 2023-11-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2052522-how-to-find-duration-of-peaks-valleys-in-chart

评论： Mathieu NOE 2023-12-11

NV1-9999.csv

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 个评论
显示无隐藏无

Dyuman Joshi 2023-11-27

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

Peter Perkins 2023-11-27

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 2023-11-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2052522-how-to-find-duration-of-peaks-valleys-in-chart#answer_1360597

编辑：Star Strider 2023-11-27

在 MATLAB Online 中打开

NV1-9999.csv

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.

.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Mathieu NOE 2023-11-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2052522-how-to-find-duration-of-peaks-valleys-in-chart#answer_1360602

在 MATLAB Online 中打开

peakseek.m

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 个评论
显示 1更早的评论隐藏 1更早的评论

Mathieu NOE 2023-11-27

在 MATLAB Online 中打开

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 2023-12-11

hello again @Sherzaad Dinah

problem solved ?

请先登录，再进行评论。

How to find duration of peaks/valleys in chart

2 个评论
显示无隐藏无

回答（2 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to find duration of peaks/valleys in chart

2 个评论 显示 无隐藏 无

回答（2 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论