how to calculate duration of peaks

18 views (last 30 days)
I have an envelope signal attached y.mat and time array t.mat
I want to calculate the duration in second between the peaks (mentionned in red color in the figure below)
How I can reach this task

Accepted Answer

Star Strider
Star Strider on 9 Jul 2022
Edited: Star Strider on 9 Jul 2022
Try this —
LD1 = load('t.mat');
ty = LD1.ty.';
LD2 = load('y.mat');
sh1 = LD2.sh1;
y = sh1.';
[pks,plocs] = findpeaks(y, 'MinPeakProminence',0.5); % Desired Peak & Index
[vys,vlocs] = findpeaks(-y, 'MinPeakProminence',0.0000001); % Valleys & Indices
for k = 1:numel(pks)
vlcs(1,k) = max(vlocs(vlocs<=plocs(k))); % First Valley Index (Beginning Of Spike)
vlcs(2,k) = min(vlocs(vlocs>=plocs(k))); % Last Valley Index (End Of Spike)
end
deltat = -diff(ty(vlcs));
figure
plot(ty,sh1)
hold on
plot(ty(plocs), pks, '^r')
plot(ty(vlcs), sh1(vlcs),'.r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
text(mean(ty(vlcs)),zeros(size(plocs)), compose('\\leftarrow \\Delta t = %.6f',deltat), 'Horiz','left', 'Vert','middle', 'Rotation',60)
producing —
EDIT — (9 Jul 2022 at 17:24)
Minor edit to my original code to make it more robust and more efficient. It produces the same reault. I have no idea how to make it more robust with respect to other data sets, since I have only one data set to work with here, and it works correctly. It will be necessary to experiment with it to get it to work with other data sets. One possibiity is to experiment with the name-value pair arguments in the second findpeaks call to be certain that it detects the beginning and end ‘valleys’.for each peak, and rejects any noise that may be in the data. My code requires that each peak has valleys on either side of it, since they are in the posted data.
EDIT — (9 Jul 2022 at 22:58)
This slightly revised code works for both sets of files —
% LD1 = load('Abdelhakim Souidi t.mat');
LD1 = load('Abdelhakim Souidi (2) t.mat');
ty = LD1.ty.';
% LD2 = load('Abdelhakim Souidi y.mat');
LD2 = load('Abdelhakim Souidi(2) y.mat');
sh1 = LD2.sh1;
y = sh1.';
[pks,plocs] = findpeaks(y, 'MinPeakProminence',0.1, 'MinPeakDistance',250); % Desired Peak & Index
[vys,vlocs] = findpeaks(-y, 'MinPeakProminence',0.0001, 'MinPeakDistance',500); % Valleys & Indices
for k = 1:numel(pks)
vlcs(1,k) = max(vlocs(vlocs<=plocs(k))); % First Valley Index
vlcs(2,k) = min(vlocs(vlocs>=plocs(k))); % Last Valley Index
end
deltat = diff(ty(vlcs));
figure
plot(ty,sh1)
hold on
plot(ty(plocs), pks, '^r')
plot(ty(vlcs), sh1(vlcs),'.r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
text(mean(ty(vlcs)),zeros(size(plocs)), compose('\\leftarrow \\Delta t = %.6f',deltat), 'Horiz','left', 'Vert','middle', 'Rotation',60)
producing —
It matters not which peak of the double-peak envelopes are identified as the actual peak, since the baseline duration is the desired result, and that simply depends on the relation of individual ‘plocs’ indices and the nearest ‘vlocs’ indices.
.
  2 Comments
Lam Ha
Lam Ha on 11 Jun 2023
I follow your code and apply it on my data but I got the problems like it. It marks the higest point with x=y. I'm looking forward to hear suggestion from u.
This is my data
Star Strider
Star Strider on 11 Jun 2023
Edited: Star Strider on 11 Jun 2023
My code does exactly what it is supposed to do with your data.
Please post this as a new Question, specifically with a description of what you want to do, and include the URL of this thread.
I will not respond to it further here, since it is not directly relevant to this thread.

Sign in to comment.

More Answers (2)

Image Analyst
Image Analyst on 9 Jul 2022
Try this:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
st = load('t.mat')
t = st.ty;
sy = load('y.mat')
sh1 = sy.sh1;
y = sh1;
plot(t, y, 'b-', 'LineWidth', 2);
grid on;
hold on;
xlabel('t', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
[peakValues, indexesOfPeaks] = findpeaks(y, 'MinPeakProminence', 0.25); % Desired Peak & Index
% Plot peaks
plot(t(indexesOfPeaks), peakValues, 'rv')
numPeaks = numel(peakValues)
% Fall down the left sides
leftIndexes = zeros(numPeaks, 1);
tLeft = zeros(numPeaks, 1);
for k = 1 : numPeaks
for k2 = indexesOfPeaks(k) : -1 : 1
thisValue = y(k2);
if thisValue < 0.1 && thisValue > y(k2+1)
leftIndexes(k) = k2;
tLeft(k) = t(leftIndexes(k));
% Put a green line at the start of the peak.
xline(tLeft(k), 'Color', 'g', 'LineWidth', 2);
break;
end
end
end
% Fall down the right sides
rightIndexes = zeros(numPeaks, 1);
tRight = zeros(numPeaks, 1);
for k = 1 : numPeaks
for k2 = indexesOfPeaks(k) : length(y)
thisValue = y(k2);
if thisValue < 0.1 && thisValue > y(k2-1)
rightIndexes(k) = k2;
tRight(k) = t(rightIndexes(k));
% Put a red line at the end of the peak.
xline(tRight(k), 'Color', 'r', 'LineWidth', 2);
break;
end
end
end
% Get the widths
peakWidths = tRight - tLeft
The left side is tLeft and is indicated by the green line.
The right side is tRight and is indicated by the red line.
The widths are given by peakWidths.
peakWidths =
0.146280579131303
0.157264103844234
0.155516724912631
0.159510733899151
0.170993509735397
0.159011482775836
0.177983025461807
0.179980029955067
0.150773839241138
0.159760359460809
  2 Comments
Abdelhakim Souidi
Abdelhakim Souidi on 11 Jul 2022
@Image Analyst you are the best. If you don't mind please how can I implement findchgpts to detect the boundaries
thank you
Image Analyst
Image Analyst on 11 Jul 2022
@Abdelhakim Souidi I don't use it that much. All I'd do is look at the documentation and the examples there and try to apply it to your data, something you can do as well as me.

Sign in to comment.


Image Analyst
Image Analyst on 9 Jul 2022
Edited: Image Analyst on 9 Jul 2022
  7 Comments
Image Analyst
Image Analyst on 11 Jun 2023
Try adjusting the 'MinPeakDistance' to make sure your valleys are not too close together.

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!