How can I calculate a mean value of peak amplitude and area under the curve?

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

[Dyuman Joshi]

Do you have to find the amplitude and AUC for every peak? or just corresponding to the global maximum peak?

[Erik Näslund]

I need to find amplitude/AUC for every peak, surrounding a given time point, and calculate a "local" mean.

Answer 1

[Star Strider]

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)

L = 449359

Fs = 960;
Fn = Fs/2;
t = linspace(0, L-1, L)/Fs;                                                                             % Time Vector
t = t(:);                                                                                               % Force Column Vector
NFFT = 2^nextpow2(L)

NFFT = 524288

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'})

PPG_Statistics = 533×5 table
    Peak    t_start    t_end     PeakToPeak       AUC    
    ____    _______    ______    __________    __________

      1     0.35625    1.0615    2.1852e+05         68335
      2      1.0615    1.9187    2.0586e+05         74229
      3      1.9187     2.701    2.8195e+05         94319
      4       2.701    3.4365    3.6124e+05    1.0798e+05
      5      3.4365    4.2542    4.7608e+05    1.4338e+05
      6      4.2542    5.1667    5.3515e+05    1.7238e+05
      7      5.1667     6.025    5.2727e+05    1.6454e+05
      8       6.025    6.7417    3.9724e+05    1.3015e+05
      9      6.7417    7.7635    3.4405e+05    1.4125e+05
     10      7.7635    8.6667    3.2516e+05    1.2902e+05
     11      8.6667     9.601    2.6599e+05    1.0231e+05
     12       9.601    10.455    2.0859e+05         75130
     13      10.455    11.417    2.6815e+05    1.0332e+05
     14      11.417    12.381     3.687e+05    1.4273e+05
     15      12.381    13.354    4.8558e+05    1.8557e+05
     16      13.354    14.235    4.9775e+05    1.7913e+05

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.

.

[Erik Näslund]

Amazing! Thank you very much. It looks like I am getting very close to what I was searching for. The summarizing table is a bonus feature. The solution is good, but If you could help me adapt it a little bit further, it would help a lot!

As an example, I am providing two MATLAB-files, consisting of two separate column vectors:

- IRn: The raw/normal IR-recording.

- IRac_filt: A bandpass filtered version of the IR-recording using a second order Butterwurth filter, with cut-off frequencys 0,4 and 30 Hz respectively.

fs = 100; %Sampling frequency in Hz
n_order = 2; % filter order
fcutlow = 0.4;   %low cut frequency in Hz
fcuthigh = 30;   %high cut frequency in Hz

[B,A] = butter(n_order,[fcutlow,fcuthigh]/(fs/2),'bandpass');

IRac_filt = filtfilt(B,A,IRn);

This particular recording is 1770 seconds long. In this example, I would like to analyse the obtained curve at six different time points (95 s, 405 s, 690 s, 985 s, 1340 s and 1625 s, scale of X-axis is 10^4, depending on sampling rate). The number of time points may vary for different recordings. At each time point, I would like to have a mean estimate of peak amplitude and AUC for instance for 6-8 seconds or 10 peaks? I don't know which method is most suitable for me?

Is there an easy fix for this?

[Star Strider]

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)

L = 176950

Fs = 100;
Fn = Fs/2;
t = linspace(0, L-1, L)/Fs;                                                                             % Time Vector
t = t(:);                                                                                               % Force Column Vector
NFFT = 2^nextpow2(L)

NFFT = 262144

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'})

PPG_Statistics = 3046×5 table
    Peak    t_start    t_end    PeakToPeak     AUC  
    ____    _______    _____    __________    ______

      1      0.46      1.02       5144.9      1497.4
      2      1.02      1.59       4546.3      1316.7
      3      1.59      2.15         4272      1209.5
      4      2.15      2.71       4409.2      1304.4
      5      2.71      3.28       5183.4      1538.4
      6      3.28      3.85       4693.4      1364.8
      7      3.85      4.41       4264.2      1217.7
      8      4.41      4.97       4384.7      1309.8
      9      4.97      5.54       5032.1      1479.1
     10      5.54      6.11       4601.4      1316.5
     11      6.11      6.67       4194.7      1198.5
     12      6.67      7.24         4356      1266.3
     13      7.24       7.8       4932.7      1454.3
     14       7.8      8.38       4629.5      1343.7
     15      8.38      8.95       4180.5      1217.4
     16      8.95      9.52       4379.6      1298.9

% 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'})

Select_PPG_Summary = 6×4 table
    t_start    t_end     Mean_P2P    Mean_AUC
    _______    ______    ________    ________

     95.41     101.13       4444      1272.1 
    405.23     410.63     4355.3      1270.8 
    690.38      695.8     4702.4      1500.6 
    985.47     990.93     5479.7      1765.1 
    1340.3     1347.4     4395.2      1470.4 
    1625.1     1630.8       3927        1254 

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]

Thank you so much! I really appreciate it. I will try this.

/Erik

[Star Strider]

As always, my pleasure!

You mentioned that you are a healthcare professional. I’m a retired B/C Internal Medicine physician.

.

[Erik Näslund]

I'm a physician specialized in anesthesia. Interesting, I guess it's never to late to learn to code :)

[Star Strider]

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]

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]

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]

Ah - It was that simple. Adjusted the factor multiplying MinPeakProminence - value to 0.05, and after that MATLAB has not, so far, misinterpreted the curve.

I have been very happy with your code so far! With this modification, it's an easy fix re-running the script for those recordings needing adjustment.

[Star Strider]

Thank you!

Answer 2

[Aman]

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!

[Erik Näslund]

Thank you very much for your answer, which is elegant. To clarify, I don't have a function to work with (?), but rather an imported column vector to MATLAB. In this vector, I need to, for X number of time points, find a mean estimate of AUC/peak amplitude at each time point.

Is your solution applicable to such a problem?

Thank you!