Identifying PTT of ECG and PPG signals

91 views (last 30 days)
Paul
Paul on 11 Apr 2024 at 15:27
Commented: Star Strider on 19 Apr 2024 at 14:57
I'm trying to detect the pulse transit time between ecg and ppg signals. I have the ecg and ppg signals filtered and normalised around 1. Here is my code so far
clc;
clear;
close all;
All_Data = {}; % set up cell array for all workbook data
Plotted_Data_ECG_PPG = {}; % Set up cell array for ECG and PPG plotted data
% Parameters for moving average filter
windowSize = 2;
ecgThreshold = 0.375;
ppgThreshold = 0.1;
for i = 1:1 % set up loop to run for 11 iterations (no. of workbooks)
filename = sprintf('%03d.xlsx', i); % set up each filename 001-011 to be read on each iteration
opt = detectImportOptions(filename); % import data from filename
sheets = sheetnames(filename); % retrieve names of sheets
for j = 1:numel(sheets) % iterate over each sheet within a workbook
tBC = readtable(filename, opt, 'Sheet', sheets{j}); % extract sheets from workbooks
Time_Data = tBC{:, 'Var4'}; % extract time data from column 4
% Remove invalid rows from table (corresponding to NaN's in Time_Data
rows2delete = isnan(Time_Data);
tBC(rows2delete,:) = [];
Time_Data = tBC{:, 'Var4'}; % extract time data from column 4
% Fit a polynomial p of 1st degree to the first points of each staircase (top curve)
xx = (1:numel(Time_Data))';
dTime_Data = diff(Time_Data);
ind = 1+find(dTime_Data>max(dTime_Data)/100);
ind = ind(2:end); % in some cases the first point is too much off trend and we prefer not to use it.
p = polyfit(xx(ind),Time_Data(ind),1);
% Evaluate the fitted polynomial p and plot:
Time_Data_new = polyval(p,xx);
% % debug plots (to double check)
% figure(j)
% plot(xx,Time_Data,xx,Time_Data_new);
ECG_Data = tBC{:, 'Var6'}; % extract ECG data from column 6
PPG_Data = tBC{:, 'Var8'}; % extract PPG data from column 8
SYS_Data = tBC{:, 'Var10'}; % extract reference systolic pressure data from column 10
DIA_Data = tBC{:, 'Var12'}; % extract reference diastolic pressure data from column 12
% Filter ECG and PPG data
ECG_Data_smoothed = movmean(ECG_Data, windowSize);
PPG_Data_smoothed = movmean(PPG_Data, windowSize);
% Normalize ECG and PPG data around 0
ECG_Data_normalized = ECG_Data_smoothed - mean(ECG_Data_smoothed);
PPG_Data_normalized = PPG_Data_smoothed - mean(PPG_Data_smoothed);
% Fnd peaks for ECG and PPG data above thresholds
[ecgPeaks, ecgPeakIdx] = findpeaks(ECG_Data_normalized, 'MinPeakHeight', ecgThreshold);
[ppgPeaks, ppgPeakIdx] = findpeaks(PPG_Data_normalized, 'MinPeakHeight', ppgThreshold);
Time_Formatted = datestr(Time_Data_new, 'HH:MM:SS:fff'); % Format time data with milliseconds
t = datetime(Time_Formatted,'InputFormat','HH:mm:ss:SSS');
t.Format = 'HH:mm:ss.SSS';
Date_Number = datenum(t);
t_plot = datetime(Date_Number, 'ConvertFrom', 'datenum');
t_plot.Format = 'HH:mm:ss.SSS';
%Plotting ECG and PPG Data on the same graph
figure;
plot(t_plot,ECG_Data_normalized, 'b', t_plot(ecgPeakIdx), ecgPeaks, 'ro', ...
t_plot, PPG_Data_normalized, 'r', t_plot(ppgPeakIdx), ppgPeaks, 'go');
xlabel('time');
ylabel('Normalized Data');
title(['ECG and PPG Data vs Time - Workbook ', num2str(i), 'sheet', sheets{j}]);
legend('ECG', 'ecgPeaks', 'PPG', 'PPG Peaks');
% Store the figure handle for later reference if needed
Plotted_Data_ECG_PPG{i,j} = gcf;
Time_Data_Cell = num2cell(Time_Data);
ECG_Data_Cell = num2cell(ECG_Data);
PPG_Data_Cell = num2cell(PPG_Data);
SYS_Data_Cell = num2cell(SYS_Data);
DIA_Data_Cell = num2cell(DIA_Data);
% Place the formatted time data and numerical data in the first row
All_Data{i,j} = [cellstr(Time_Formatted), ECG_Data_Cell, PPG_Data_Cell, SYS_Data_Cell, DIA_Data_Cell];
end
end
  1 Comment
Manikanta Aditya
Manikanta Aditya on 11 Apr 2024 at 15:54
Hi,
Your code seems to be well-structured for processing ECG and PPG signals. However, to calculate the Pulse Transit Time (PTT), you need to find the time difference between the corresponding peaks of the ECG and PPG signals.

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Star Strider
Star Strider on 11 Apr 2024 at 16:07
Without your data, it is going to be nearly impossible to suggest an approach to this. Also, without your data, it is impossible to understand how your code works with it.
First, there are several ways to eliminate baseline variation and offset. One is a simple highpass filter. (Use 'ImpulseResponse','iir' for best results.) However that may not always be the best option, since it can produce artifacts. Another way is to us the islocalmin function to detect minima along the trace, fit a low-order polynomial to those poiints, and then subtract a vector using that function and representing the entire signal (usually the time vector) from the signal.
Second, the R deflection is always going to preceed the succeeding PPG peak, so take the minimum time value for the PPG peak that is greater than a specific R-deflection time to detect it. The min function should be able to do this. If you know the approximate distance from the heart to the PPG sensor, you can possibly calculate a propagation velocity.
.
  7 Comments
Show 5 older comments
Kern
Kern on 19 Apr 2024 at 14:43
That does make sense, I was trying to use the resample function on the signal processing toolbox before, but I hadn't thought of using it with the time values aswell, if I am correct in understanding. Thank you anyway, I will let you know how I get on.
By PAT, pulse arrival time, I just mean PTT, I have seen it used interchangably in the literature.
Star Strider
Star Strider on 19 Apr 2024 at 14:57
My pleasure!
Thank you. I’m not familiar with all the literature, and haven’t worked with this sort of data in a while.
I usually use this syntax Resample Nonuniformly Sampled Data in Timetable for resample with a timetable (that would be the preferred approach here). It should be relatively straightforward to extract the sampling frequency afterwards, if you need that for filter design or other applications.

Sign in to comment.

More Answers (0)

Categories

Find more on Time-Frequency Analysis in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!