필터 지우기
필터 지우기

How to find duration of peaks/valleys in chart

조회 수: 9 (최근 30일)
Sherzaad Dinah
Sherzaad Dinah 2023년 11월 27일
댓글: Mathieu NOE 2023년 12월 11일
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
Dyuman Joshi 2023년 11월 27일
What does 'duration' mean in this context? Please specify.
Peter Perkins
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.

댓글을 달려면 로그인하십시오.

답변 (2개)

Star Strider
Star Strider 2023년 11월 27일
편집: Star Strider 2023년 11월 27일
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.
.

Mathieu NOE
Mathieu NOE 2023년 11월 27일
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
Mathieu NOE
Mathieu NOE 2023년 11월 27일
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
Mathieu NOE 2023년 12월 11일
hello again @Sherzaad Dinah
problem solved ?

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Multirate Signal Processing에 대해 자세히 알아보기

제품


릴리스

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by