How to calculate 50% and 66% width of each peak in the data?

Hi

I have the PPG data collected at 1kHz sampling freuqncy. I need to store the 33%, 50% and 66% width time duration(x-axis values, which in total is of 2.1s) of each detected peaks.

Please help me how i can do that.?

Answer 1

[Star Strider]

Try this —

LD = load('matlab.mat');
data = LD.filteredData.';

L = numel(data)

L = 2100

Fs = 1000;
x = linspace(0, L-1, L).'/Fs;

figure
plot(x, data)
grid
title('Original Data')

[pks,plocs] = findpeaks(data, 'MinPeakProminence',0.25);
[vys,vlocs] = findpeaks(-data, 'MinPeakProminence',0.25);
vlocs = [1; vlocs; numel(x)];
vys = data(vlocs);

amplvls = [1/3; 1/2; 2/3];
for k = 1:numel(pks)
    [~,idx] = mink(abs(plocs(k)-vlocs),2);                                      % Indices Of Two Nearest Valleys To Peak
    idxv(:,k) = sort(vlocs(idx));                                               % Sort Them
    idxrng = idxv(1,k) : idxv(2,k);                                             % Index Range
    midpt(k) = x(plocs(k));                                                     % 'x' Value For Peak
    ymp(k) = [midpt(k), 1] * ([x(idxv(:,k)) ones(2,1)] \ data(idxv(:,k)));      % 'y' Value Of Interpolated Baseline For This Peak
    pkamp(k) = pks(k) - ymp(k);                                                 % Peak Value With Respect To Baseline
    idxrng1 = idxv(1,k) : plocs(k);                                             % Index Range: First Valley To Peak
    idxrng2 = plocs(k) : idxv(2,k);                                             % Index Range: Peak To Second Valley
    t_start(:,k) = interp1(data(idxrng1), x(idxrng1), pkamp(k)*amplvls);        % Start Tiimes
    t_end(:,k) = interp1(data(idxrng2), x(idxrng2), pkamp(k)*amplvls);          % End Times
    t_dif(:,k) = t_end(:,k) - t_start(:,k);                                     % Time Difference
end

LevelTimes = array2table(t_dif, 'VariableNames',compose('Peak %2d',1:numel(pks)), 'RowNames',compose('%.0f%%',amplvls*100))

LevelTimes = 3×3 table
           Peak  1    Peak  2    Peak  3
           _______    _______    _______

    33%    0.23376    0.22093    0.23351
    50%    0.18972     0.1792    0.19022
    67%    0.15338    0.13833    0.14962

figure
hp{1} = plot(x, data, 'DisplayName','Data');
hold on
for k = 1:numel(pkamp)
    hp{2} = plot(x(idxv(:,k)), data(idxv(:,k)), '--k', 'DisplayName','Baselines');
    hp{3} = plot([1 1]*midpt(k), [ymp(k) pks(k)], '-r', 'DisplayName','Peak Amplitudes Above Inter-Peak Baseline');
    hp{4} = plot([t_start(:,k) t_end(:,k)].', (pkamp(k)*amplvls*[1 1]).', '-g', 'LineWidth',1.5);
end
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude')
title({'Original Data With Peak Heights Above Baselilne','and Level Lines'})
% legend([hp{1} hp{2}(1) hp{3}(1)], 'Location','best')
text(midpt, ymp+pkamp/2, compose('\\leftarrow %.2f',pkamp))
for k = 1:numel(pks)
text(t_start(:,k), amplvls*pkamp(k), compose('%.3f s \\rightarrow',t_dif(:,k)), 'Horiz','right', 'Vert','middle', 'FontSize',9)
end

With more PPG pulses, the time annotations will not be able to be displayed simply because there will not be room for all of them. They work here.

.

[Star Strider]

The first approach uses the amplitude of the peak above the baseline to define the levels.

Defining different fractional levels for the rising and falling limbs of the PPG pulse —

LD = load('matlab.mat');
data = LD.filteredData.';

L = numel(data);
Fs = 1000;
x = linspace(0, L-1, L).'/Fs;

figure
plot(x, data)
grid
title('Original Data')

[pks,plocs] = findpeaks(data, 'MinPeakProminence',0.25);
[vys,vlocs] = findpeaks(-data, 'MinPeakProminence',0.25);
vlocs = [1; vlocs; numel(x)];
vys = data(vlocs);

amplvls = [1/3; 1/2; 2/3];
for k = 1:numel(pks)
    [~,idx] = mink(abs(plocs(k)-vlocs),2);                                      % Indices Of Two Nearest Valleys To Peak
    idxv(:,k) = sort(vlocs(idx));                                               % Sort Them
    idxrng = idxv(1,k) : idxv(2,k);                                             % Index Range
    midpt(k) = x(plocs(k));                                                     % 'x' Value For Peak
    ymp(k) = [midpt(k), 1] * ([x(idxv(:,k)) ones(2,1)] \ data(idxv(:,k)));      % 'y' Value Of Interpolated Baseline For This Peak
    pkamp(k) = pks(k) - ymp(k);                                                 % Peak Value With Respect To Baseline
    idxrng1 = idxv(1,k) : plocs(k);                                             % Index Range: First Valley To Peak
    idxrng2 = plocs(k) : idxv(2,k);                                             % Index Range: Peak To Second Valley
    amp1(:,k) = data(plocs(k)) - data(idxv(1,k));                               % Rise Amplitude
    amp2(:,k) = data(plocs(k)) - data(idxv(2,k));                               % Fall Amplitude
    ampmtx(:,:,k) = amplvls*[amp1(:,k) amp2(:,k)]+data(idxv(:,k)).';            % Matrix Of Fractional Amplitudes For Each Secton Of Each PPG Pulse
    t_start(:,k) = interp1(data(idxrng1), x(idxrng1), ampmtx(:,1,k));           % Start Tiimes
    t_end(:,k) = interp1(data(idxrng2), x(idxrng2), ampmtx(:,2,k));             % End Times
    t_dif(:,k) = t_end(:,k) - t_start(:,k);                                     % Time Difference
end

LevelTimes = array2table(t_dif, 'VariableNames',compose('Peak %2d',1:numel(pks)), 'RowNames',compose('%.0f%%',amplvls*100))

LevelTimes = 3×3 table
           Peak  1    Peak  2    Peak  3
           _______    _______    _______

    33%    0.21272    0.21444    0.21964
    50%    0.17348    0.17316    0.17807
    67%    0.13441    0.13256    0.13679

figure
hp{1} = plot(x, data, 'DisplayName','Data');
hold on
for k = 1:numel(pkamp)
    hp{2} = plot(x(idxv(:,k)), data(idxv(:,k)), '--k', 'DisplayName','Baselines');
    hp{3} = plot([1 1]*midpt(k), [ymp(k) pks(k)], '-r', 'DisplayName','Peak Amplitudes Above Inter-Peak Baseline');
    hp{4} = plot([t_start(:,k) t_end(:,k)].', [ampmtx(:,1,k) ampmtx(:,2,k)].', '-g', 'LineWidth',1.5);
end
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude')
title({'Original Data With Peak Heights Above Baselilne','and Level Lines'})
% legend([hp{1} hp{2}(1) hp{3}(1)], 'Location','best')
text(midpt, ymp+pkamp*0.75, compose('\\leftarrow %.2f',pkamp))
for k = 1:numel(pks)
text(t_start(:,k), ampmtx(:,1,k), compose('%.3f s \\rightarrow',t_dif(:,k)), 'Horiz','right', 'Vert','middle', 'FontSize',9)
end

Since the rising and falling limbs have different amplitudes, this version calculates the fractional amplitudes for each limb, and the corresponding times and elapsed times.

.