How can I calculate a mean value of peak amplitude and area under the curve?
Show older comments
Hello,
I hope someone can help me with my problem. As a healthcare professional, I lack the basic understanding of Data programming and mathematics.
I have a recording of, how the absorption varies, from an infrared photoplethysmography signal obtained in vivo. The data is imported in Matlab as column vectors. I have managed to apply a low pass and a bandpass filter to the raw signal, separating the different components I have to analyse further. From the band pass filtered signal (picture below), I need to – for a number of given time points (illustrated with a green arrow):
- Determine the peak amplitude (local minima to the peak).
- Determine the area under the curve (AUC) – from minima to the peak.
- The amplitude and AUC should be estimated with some mean value, e.g. +/- 5 seconds from the given time point or +/- 5 peaks. I don’t know which of them is the best method?
Does anyone have any suggestions for possible solutions?
I appreciate any help you can provide.
Y-axis: light absorption in arbitrary units, X-axis: time, based on sampling rate (100 Hz) – every second = 100 units
2 Comments
Dyuman Joshi
on 15 Feb 2023
Do you have to find the amplitude and AUC for every peak? or just corresponding to the global maximum peak?
Erik Näslund
on 15 Feb 2023
Accepted Answer
It always helps to have a signal to work with, so I found one from a previous pulse-plethysmograph post and am using it here. (All the pre-processing steps may not be necessary for your signal, however they were for this one. They also demonstrate how to pre-process a physiological signal in general.)
Try something like this —
PPG = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/103608/p3_sit.csv'); % Use PPG Signal From Previoous Answer
PPG = PPG(PPG>=1E+6); % Remove Artefacts (This Signal Only)
L = size(PPG,1)
Fs = 960;
Fn = Fs/2;
t = linspace(0, L-1, L)/Fs; % Time Vector
t = t(:); % Force Column Vector
NFFT = 2^nextpow2(L)
FT_PPG = fft((PPG-mean(PPG)).*hann(L),NFFT);
Fv = linspace(0,1,NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FT_PPG(Iv)*2))
grid
xlim([0 10])
PPG = lowpass(PPG, 5, Fs, 'ImpulseResponse','iir'); % Remove High-Frequency Noise (This Signal Only)
PPGmax = max(PPG);
PPGmin = min(PPG);
[pks,plocs] = findpeaks(PPG, 'MinPeakProminence',(PPGmax-PPGmin)*0.05); % Peaks
[vys,vlocs] = findpeaks(-PPG, 'MinPeakProminence',(PPGmax-PPGmin)*0.05); % Valleys
for k = 1:numel(vlocs)-1
idxrng = vlocs(k) : vlocs(k+1); % Index Range For This Pulse
idxlim = [vlocs(k); vlocs(k+1)]; % Inmdex Limits
tv = t(idxrng); % Pulse 't' Vector
sv = PPG(idxrng); % Pulse 'PPG' Vector
ts(k,:) = t(vlocs(k)); % Pulse Start Time
te(k,:) = t(vlocs(k+1)); % Pul;se End Time
kv(k,:) = k; % Peak Counter
B = [t(idxlim) ones(2,1)] \ PPG(idxlim); % Linear Detrending Regression Parameters
dtsv = [t(idxrng) ones(size(t(idxrng)))] * B; % Linear Detrending Vector
svc = sv - dtsv; % Detrended Signal
p2p(k,:) = max(svc); % Peak-To-Peak Amplitude
AUC(k,:) = trapz(tv,svc); % Pulse Area
end
PPG_Statistics = table(kv,ts,te,p2p,AUC, 'VariableNames',{'Peak','t_start','t_end','PeakToPeak','AUC'})
figure
hp{1} = plot(t, PPG, 'DisplayName','PPG Data');
hold on
hp{2} = plot(t(plocs), pks, '^r', 'DisplayName','Peaks');
hp{3} = plot(t(vlocs), -vys, 'vr', 'DisplayName','Valleys');
hp{4} = plot((t(plocs(1:numel(p2p)))*[1 1]).', [pks(1:numel(p2p)) pks(1:numel(p2p))-p2p].', '-k', 'DisplayName','Peak-Peak');
hp{5} = plot(t(vlocs), -vys, '-g', 'DisplayName','Inter-Valley Baseline');
hold off
grid
legend([hp{1}(1),hp{2}(1),hp{3}(1),hp{4}(1),hp{5}(1)],'Location','best')
xlim([0 10]) % Zoom To See Details
I have no idea what the original units were for this signal, however the values for the derived statistics appear to be correct.
Try my code with your signal. If there are any problems, I will help adapt my code to it.
.
10 Comments
Erik Näslund
on 15 Feb 2023
I used the unfiltered data, since my code is relatively robust to baseline variations and high-frequency noise.
With respect to the data at those points, I created the ‘tiv’ vector to guide the subsequent analysis.
Getting the information on the mean values over several consecutive PPG pulse spans from the ‘PPG_Statistics’ table is straightforward. The ‘Select_PPG_Summary’ table uses the ‘PPG_Statistics’ table to initially find the first ‘t_start’ value meeting the criteria in ‘tiv’, and then calculates the mean of the ‘PeakToPeak’ and ‘AUC’ values for that peak and the next nine rows, for a total to 10 consecutive peaks.
Try this —
% PPG = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/103608/p3_sit.csv'); % Use PPG Signal From Previoous Answer
LD = load(websave('IRn','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1296730/IRn.mat'));
PPG = LD.IRn; % Remove Artefacts (This Signal Only)
L = size(PPG,1)
Fs = 100;
Fn = Fs/2;
t = linspace(0, L-1, L)/Fs; % Time Vector
t = t(:); % Force Column Vector
NFFT = 2^nextpow2(L)
FT_PPG = fft((PPG-mean(PPG)).*hann(L),NFFT);
Fv = linspace(0,1,NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FT_PPG(Iv)*2))
grid
xlim([0 10])
xlabel('Frequency (Hz)')
ylabel('Magnitude (Units)')
PPG = lowpass(PPG, 10, Fs, 'ImpulseResponse','iir'); % Remove High-Frequency Noise (This Signal Only)
PPGmax = max(PPG);
PPGmin = min(PPG);
[pks,plocs] = findpeaks(PPG, 'MinPeakProminence',(PPGmax-PPGmin)*0.01); % Peaks
[vys,vlocs] = findpeaks(-PPG, 'MinPeakProminence',(PPGmax-PPGmin)*0.01); % Valleys
if plocs(1) < vlocs(1)
plocs = plocs(2:end);
pks = pks(2:end);
end
% plocs
% vlocs
for k = 1:numel(vlocs)-1
idxrng = vlocs(k) : vlocs(k+1); % Index Range For This Pulse
idxlim = [vlocs(k); vlocs(k+1)]; % Inmdex Limits
tv = t(idxrng); % Pulse 't' Vector
sv = PPG(idxrng); % Pulse 'PPG' Vector
ts(k,:) = t(vlocs(k)); % Pulse Start Time
te(k,:) = t(vlocs(k+1)); % Pul;se End Time
kv(k,:) = k; % Peak Counter
B = [t(idxlim) ones(2,1)] \ PPG(idxlim); % Linear Detrending Regression Parameters
dtsv = [t(idxrng) ones(size(t(idxrng)))] * B; % Linear Detrending Vector
svc = sv - dtsv; % Detrended Signal
p2p(k,:) = max(svc); % Peak-To-Peak Amplitude
AUC(k,:) = trapz(tv,svc); % Pulse Area
end
PPG_Statistics = table(kv,ts,te,p2p,AUC, 'VariableNames',{'Peak','t_start','t_end','PeakToPeak','AUC'})
% timespan = [min(ts) max(te)]
figure
hp{1} = plot(t, PPG, 'DisplayName','PPG Data');
hold on
hp{2} = plot(t(plocs), pks, '^r', 'DisplayName','Peaks');
hp{3} = plot(t(vlocs), -vys, 'vr', 'DisplayName','Valleys');
hp{4} = plot((t(plocs(1:numel(p2p)))*[1 1]).', [pks(1:numel(p2p)) pks(1:numel(p2p))-p2p].', '-k', 'DisplayName','Peak-Peak');
hp{5} = plot(t(vlocs), -vys, '-g', 'DisplayName','Inter-Valley Baseline');
hold off
grid
legend([hp{1}(1),hp{2}(1),hp{3}(1),hp{4}(1),hp{5}(1)],'Location','best')
xlim([0 10]) % Zoom To See Details
xlabel('Time (s)')
ylabel('Amplitude (Units)')
tiv = [95; 405; 690; 985; 1340; 1625];
% PPGiv = interp1(t, PPG, tiv)
for k = 1:numel(tiv) % Starting Time
Rows = find(PPG_Statistics.t_start >= tiv(k),1) + (0:9); % Next 10 Rows
Rows = Rows(1):min(Rows(end),size(PPG_Statistics,1)); % Truncate To End Of 'PPG_STtatistics' If Necessary
tst(k,:) = PPG_Statistics.t_start(Rows(1));
tnd(k,:) = PPG_Statistics.t_end(Rows(end)); % Ending Time
Result(k,:) = mean(PPG_Statistics{Rows,[4 5]}); % Calculate Mean Values
end
Select_PPG_Summary = table(tst,tnd,Result(:,1),Result(:,2), 'VariableNames',{'t_start','t_end','Mean_P2P','Mean_AUC'})
This should be reasonably robust for similar data, since it does a reasonable amount of checking to be certain that it does not exceed the limits of the data or the derived table arrays. Fortunately, PPG data is relatively straightforward to work with, providing that it is not excessively noisy, especially with broadband noise. So long as the noise is band-limited, it should be relatively easy to work with the signals. .
.
Erik Näslund
on 15 Feb 2023
Star Strider
on 15 Feb 2023
As always, my pleasure!
You mentioned that you are a healthcare professional. I’m a retired B/C Internal Medicine physician.
.
Erik Näslund
on 16 Feb 2023
Star Strider
on 16 Feb 2023
I seriously considered Anesthesiology, since I did an elective in it in medical school and liked it. I preferred Internal Medicine, and after my residency did a fellowship in Endocrinology, although never bothered to sit for my boards. (I would have preferred Intensive Care medicine to Endocrinology, however it was only beginning to be available then.) Instead, I went back to graduate school and earned a M.S. in Biomedical Engineering, since personal computers were becoming available and I wanted to learn how to make the most of that technology in a medical context. (That was one of my best life decisions.) My subsequent career ended up combining both disciplines.
Erik Näslund
on 14 Jun 2023
Mr Star Strider,
Again, thank you for your solution it has help me a lot. However, from time to time MATLAB misinterprets an artifact or a notch for a peak, resulting in a value less than expected, since the values presented is a mean of 10 consecutive peaks. I have managed to recalculate P2P, simply by manually selecting 3 consecutive peaks and calculating a new mean, using the values from the peak and valleys. However, I have not been equally successfull in recalculating AUC. I have been using the trapz-function but have not been able to subtract the area between the Inter-Valley baseline and y-axis = 0.
I would be very happy if you could help me with this!
Sincerely
/Erik
In the example picture below - I would like to find the AUC for the peaks between 128,32 s - 130,5 s.
Star Strider
on 14 Jun 2023
You will need to adjust the 'MinPeakProminence' value (probably increase them) in at least the first (and perhaps both) findpeaks calls.
I did my best to make my code generally applicable, however one size does not fit all!
Erik Näslund
on 14 Jun 2023
Star Strider
on 14 Jun 2023
Thank you!
More Answers (1)
Aman
on 15 Feb 2023
Hi,
I understand that you have a function, and you want to find out its global maxima. Then you want to find the nearest local minima on the either side of the global maxima and then calculate the area between the two local minima.
For finding the global maxima you can use the “max” function. Upon getting the global maxima you can use the “islocal” function to get all the local minima, upon getting all the local minima you can find the nearest local minima on the either side of the global maxima. Once you get the local minima of the either side you can find the subsection between them and then use the ”trapz” function to get the area between them. You can refer to the following example:
x = 1:100;
A = (1-cos(2*pi*0.01*x)).*sin(2*pi*0.15*x); % objective function
[maxVal,indx] = max(A); % global maxima
TF = islocalmin(A); % all local minima
for i=indx:-1:1
if(TF(i))
localMinLeft = i;
localMinLeftVal = A(i); % local minima to the left of the global maxima
break;
end
end
for i=indx:100
if(TF(i))
localMinRight = i;
localMinRightVal = A(i); % local minima to the right of the global maxima
break;
end
end
subSection = A(localMinLeft:localMinRight); % subsection between the local minimas
plot(x,A);
hold on;
plot(localMinRight,localMinRightVal,'r*');
plot(localMinLeft,localMinLeftVal,'r*');
plot(x(indx),A(indx),'bo');
plot(localMinLeft:localMinRight,subSection);
hold off;
disp('Area under the curve is ' + trapz(subSection)); % area under the subsection selected
Hope this helps!
1 Comment
Erik Näslund
on 15 Feb 2023
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
