%% Data Processing for 2023 ARM Trial
% Ella Guy
% Last Updated: 24FEB2024
clear all, close all, clc
%--------------------------------------------------------------------------
% Inputs ------------------------------------------------------------------
SubjectNumber = 20; % de-identifyed subject number
TestNumber = 2; % test reference (1: increasing PEEP,
% 2: Breath-holds with increasing PEEP, and
% 3: Forced Expiratory Manuevers)
breathMin = 0.007; % minimum fluctuation in second derivative breath selection
% Alignment Breath --------------------------------------------------------
% In BOB data
T_min_BOB = 26; %[s]
T_max_BOB = 27; %[s]
% In EIT data
T_min_EIT = 26; %[s]
T_max_EIT = 27; %[s]
%--------------------------------------------------------------------------
%% Loads Data
%--------------------------------------------------------------------------
Trials = {'PEEP', 'PEEP_BH', 'FEM'};
Trial = Trials{TestNumber};
RecordingLengths = [350, 350, 60];
RecordingLength = RecordingLengths(TestNumber);
% sensor offsets
OffsetA = -0.150;
OffsetB = -0.064;
%--------------------------------------------------------------------------
% Load Respiratory Data ---------------------------------------------------
Infile_format = 'PQ_rawData/Subject%d_%s.csv';
Infile = sprintf(Infile_format, SubjectNumber, Trial);
Respiratory = readtable(Infile,'NumHeaderLines',1);
time = table2array(Respiratory(:,1)); % time [s]
GaugeP = table2array(Respiratory(:,2)); % Gauge pressure reading [cmH2O]
InhaleDeltaP = table2array(Respiratory(:,3)); % Inspiratory differenrtial pressure [cmH2O]
ExhaleDeltaP = table2array(Respiratory(:,4)); % Expiratory differenrtial pressure [cmH2O]
Depth_chest = table2array(Respiratory(1,5)); % Chest crossectional depth (a) [cm]
Width_chest = table2array(Respiratory(1,6)); % Chest crossectional width (b) [cm]
StartTime = datestr(table2array(Respiratory(1,7)));
BOB_timestep = 0.01;
[h, m, s] = hms(datetime(StartTime, 'InputFormat', 'HH:mm:ss'));
PSD_StartTime = [h, m, s];
for ai = length(time):-1:1
if time(ai) >= 0
trim_PSG = ai;
break
end
end
PSD_t = time(1:trim_PSG);
PSD_P = GaugeP(1:trim_PSG);
PSD_DPI = InhaleDeltaP(1:trim_PSG);
PSD_DPE = ExhaleDeltaP(1:trim_PSG);
% -------------------------------------------------------------------------
% Load EIT File -----------------------------------------------------------
EIT_infile_format = 'EIT_rawData/S%02d_%s.bin';
EIT_infile = sprintf(EIT_infile_format, SubjectNumber, Trial);
EIT_data = read_binData(EIT_infile);
Aeration = sum(sum(EIT_data));
Aeration = Aeration(:);
EIT_sample_frequency = 50; % [Hz]
EIT_timestep = 1/EIT_sample_frequency; % [s]
Time_Aeration = EIT_timestep.*([1:length(Aeration)]);
% -------------------------------------------------------------------------
% Load ECG File -----------------------------------------------------------
ECG_infile_format = 'HRM_rawData/ECG/S%02d_ECG.csv';
ECG_infile = sprintf(ECG_infile_format, SubjectNumber);
ECG_file = readtable(ECG_infile,'NumHeaderLines',1);
ECG = table2array(ECG_file(:,2));
ECG_time = table2array(ECG_file(:,1));
[h1, m1, s1] = hms(datetime(ECG_time(1), 'InputFormat', 'HH:mm:ss.SSS'));
ECG_StartTime = [h1, m1, s1];
trim_ECG = 1;
for ei = 1:length(ECG_time)
[hw, mw, sw] = hms(datetime(ECG_time(ei), 'InputFormat', 'HH:mm:ss.SSS'));
if hw == PSD_StartTime(1)
if mw == PSD_StartTime(2)
if floor(sw) == PSD_StartTime(3)
trim_ECG = ei;
break
end
end
end
end
ECG_time2 = ECG_time(trim_ECG:end);
for bii = 1:length(ECG_time2)
[hw, mw, sw] = hms(datetime(ECG_time2(bii), 'InputFormat', 'HH:mm:ss.SSS'));
ECG_t(bii) = (hw - PSD_StartTime(1))*3600 + ...
(mw - PSD_StartTime(2))*60 + ...
(sw - PSD_StartTime(3));
if ECG_t(bii) >= RecordingLength
break
end
end
ECG_trace = ECG(trim_ECG: trim_ECG+length(ECG_t)-1);
% -------------------------------------------------------------------------
% Load PPG File -----------------------------------------------------------
PPG_infile_format = 'HRM_rawData/PPG/S%02d_PPG.csv';
PPG_infile = sprintf(PPG_infile_format, SubjectNumber);
PPG_file = readtable(PPG_infile,'NumHeaderLines',1);
PPG0 = table2array(PPG_file(:,2));
PPG1 = table2array(PPG_file(:,3));
PPG2 = table2array(PPG_file(:,4));
PPG_time = table2array(PPG_file(:,1));
[h2, m2, s2] = hms(datetime(PPG_time(1), 'InputFormat', 'HH:mm:ss.SSS'));
PPG_StartTime = [h2, m2, s2];
trim_PPG = 1;
for ci = 1:length(PPG_time)
[hw, mw, sw] = hms(datetime(PPG_time(ci), 'InputFormat', 'HH:mm:ss.SSS'));
if hw == PSD_StartTime(1)
if mw == PSD_StartTime(2)
if floor(sw) == PSD_StartTime(3)
trim_PPG = ci;
break
end
end
end
end
PPG_time2 = PPG_time(trim_PPG:end);
for cii = 1:length(PPG_time2)
[hw, mw, sw] = hms(datetime(PPG_time2(cii), 'InputFormat', 'HH:mm:ss.SSS'));
PPG_t(cii) = (hw - PSD_StartTime(1))*3600 + ...
(mw - PSD_StartTime(2))*60 + ...
(sw - PSD_StartTime(3));
if PPG_t(cii) >= RecordingLength
break
end
end
PPG_trace0 = PPG0(trim_PPG: trim_PPG+length(PPG_t)-1);
PPG_trace1 = PPG1(trim_PPG: trim_PPG+length(PPG_t)-1);
PPG_trace2 = PPG2(trim_PPG: trim_PPG+length(PPG_t)-1);
% -------------------------------------------------------------------------
% Load HR Belt File -------------------------------------------------------
HRB_infile_format = 'HRM_rawData/HRB/%d.txt';
HRB_infile = sprintf(HRB_infile_format, SubjectNumber);
HRB_file = readtable(HRB_infile,'NumHeaderLines',1);
HRB_time_str = table2array(HRB_file(:,1));
HRB_msgid = table2array(HRB_file(:,2));
HRB_rr_ms = table2array(HRB_file(:,3));
HRB_hr_bpm = table2array(HRB_file(:,4));
for i = 1:length(HRB_time_str)
HRB_time(i) = datetime(HRB_time_str{i},'InputFormat', 'HH:mm:ss.SSSSSS''Z');
end
[h3, m3, s3] = hms(HRB_time(1));
HRB_StartTime = [h3, m3, s3];
if PSD_StartTime(1) >= 12
PSD_StartTime(1) = PSD_StartTime(1) - 12;
end
trim_HRB = 1;
for di = 1:length(HRB_time)
[hw, mw, sw] = hms(HRB_time(di));
if hw >= 12
hw = hw - 12;
end
if hw == PSD_StartTime(1)
if mw == PSD_StartTime(2)
if floor(sw) == PSD_StartTime(3)
trim_HRB = di;
break
end
end
end
end
HRB_time2 = HRB_time(trim_HRB:end);
for dii = 1:length(HRB_time2)
[hw, mw, sw] = hms(HRB_time2(dii));
if hw > 12
hw = hw - 12;
end
HRB_t(dii) = (hw - PSD_StartTime(1))*3600 + ...
(mw - PSD_StartTime(2))*60 + ...
(sw - PSD_StartTime(3));
if HRB_t(dii) >= RecordingLength
break
end
end
HRB_HR = HRB_hr_bpm(trim_HRB: trim_HRB+length(HRB_t)-1);
HRB_RR = (trim_HRB: trim_HRB+length(HRB_t)-1);
% -------------------------------------------------------------------------
% Processing Respiratory Data --------------------------------------------
% Flow --------------------------------------------------------------------
d1 = 0.015; % Venturi intake diameter [m]
d2 = 0.010; % Venturi constriction diameter [m]
cd = 0.97; % Discharge coefficient
rho = 1.225; % Density of air [kg/m^3]
A1 = pi*(d1/2)^2; % Venturi intake cross-sectional area [m^2]
A2 = pi*(d2/2)^2; % Venturi constrictioncross-sectional area [m^2]
DeltaP = InhaleDeltaP(1:trim_PSG) - ExhaleDeltaP(1:trim_PSG);
PSD_Q = zeros(length(DeltaP), 1);
for ei = 1 : length(DeltaP)
if PSD_P(ei) <= 0
PSD_Q(ei) = 1000*cd*A2*sqrt((2*(abs(DeltaP(ei))*98.0665))/(rho*(1-(A2/A1)^2)));
else
PSD_Q(ei) = -1.*1000*cd*A2*sqrt((2*(abs(DeltaP(ei))*98.0665))/(rho*(1-(A2/A1)^2)));
end
end
% -------------------------------------------------------------------------
% Indicies ----------------------------------------------------------------
Smoothed_Q = movmedian(PSD_Q, 20);
Smoothed_P = smooth(movmedian(PSD_P, 20), 250);
Inflect_P = smooth(movmedian(diff(Smoothed_P,2),20), 250);
InspInd_working1 = [];
for fi = 2:length(Inflect_P)
if Inflect_P(fi-1) < 0 && Inflect_P(fi) >= 0
InspInd_working1 = [InspInd_working1, fi];
end
end
InspInd_working2 = [InspInd_working1(1)];
for fii = 2:length(InspInd_working1)-1
breath = max(cumtrapz(Inflect_P(InspInd_working1(fii):InspInd_working1(fii+1))));
if breath >= breathMin
InspInd_working2 = [InspInd_working2, InspInd_working1(fii)];
end
end
InspInd = [InspInd_working2(1)];
for fiii = 2:length(InspInd_working2)
breath = min(cumtrapz(Inflect_P(InspInd_working2(fiii):-1:InspInd_working2(fiii-1))));
if breath <= -1*breathMin
InspInd = [InspInd, InspInd_working2(fiii)];
end
end
% -------------------------------------------------------------------------
% Volume ------------------------------------------------------------------
V_working = 0;
PSD_V = [V_working];
for hi = 2:length(PSD_Q)
if ismember(hi, InspInd) == 0
V_working = V_working + trapz(BOB_timestep, PSD_Q(hi-1:hi));
else
V_working = 0;
end
PSD_V = [PSD_V ; V_working];
end
% -------------------------------------------------------------------------
figure(1)
subplot(2,1,1)
hold on
plot(PSD_t, PSD_V)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(2,1,2)
hold on
plot(Time_Aeration, Aeration)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
%% Aligning EIT Data ------------------------------------------------------
I_min_BOB = (BOB_timestep*floor(T_min_BOB/BOB_timestep))/BOB_timestep;
I_max_BOB = (BOB_timestep*floor(T_max_BOB/BOB_timestep))/BOB_timestep;
I_min_EIT = (EIT_timestep*floor(T_min_EIT/EIT_timestep))/EIT_timestep;
I_max_EIT = (EIT_timestep*floor(T_max_EIT/EIT_timestep))/EIT_timestep;
[LM_BOB, LMI_BOB] = max(PSD_V(I_min_BOB:I_max_BOB));
[LM_EIT, LMI_EIT] = max(Aeration(I_min_EIT:I_max_EIT));
LMT_BOB = PSD_t(I_min_BOB + LMI_BOB);
LMT_EIT = Time_Aeration(I_min_EIT + LMI_EIT);
% Aligning ----------------------------------------------------------------
T_align = LMT_BOB - LMT_EIT;
EIT_time = Time_Aeration + T_align;
for ji = 1:length(EIT_time)
if EIT_time(ji)>=0
trim1_EIT = ji;
break
end
end
trim2_EIT = length(EIT_time);
for jii = 1:length(EIT_time)
if EIT_time(jii)>=RecordingLength
trim2_EIT = jii;
break
end
end
EIT_t = EIT_time(trim1_EIT:trim2_EIT);
EIT_A = Aeration(trim1_EIT:trim2_EIT);
% -------------------------------------------------------------------------
% Plots ------------------------------------------------------------------
figure(2)
subplot(7,1,1)
hold on
plot(PSD_t, PSD_P)
plot(PSD_t(InspInd), PSD_P(InspInd), '.')
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(7,1,2)
hold on
plot(PSD_t, PSD_Q)
plot(PSD_t(InspInd), PSD_Q(InspInd), '.')
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(7,1,3)
hold on
plot(PSD_t, PSD_V)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(7,1,4)
hold on
plot(EIT_t, EIT_A)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(7,1,5)
hold on
plot(ECG_t, ECG_trace)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(7,1,6)
hold on
plot(PPG_t, PPG_trace0)
plot(PPG_t, PPG_trace1)
plot(PPG_t, PPG_trace2)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
subplot(7,1,7)
hold on
plot(HRB_t, HRB_HR)
xline(PSD_t(InspInd), '--')
xlim([0 RecordingLength])
grid on
grid minor
hold off
%% Saves Data
csv_outfile_format = 'ProcessedDataset/ProcessedData_Subject%02d_%s.csv';
csv_outfile_name = sprintf(csv_outfile_format, SubjectNumber, Trial);
Output = zeros(length(PPG_t), 18);
Output(1:length(PSD_t),1) = PSD_t;
Output(1:length(PSD_t),2) = PSD_P;
Output(1:length(PSD_t),3) = PSD_DPI;
Output(1:length(PSD_t),4) = PSD_DPE;
Output(1:length(PSD_t),5) = PSD_Q;
Output(1:length(PSD_t),6) = PSD_V;
Output(1:length(InspInd), 7) = InspInd;
Output(1:length(EIT_t), 8) = EIT_t;
Output(1:length(EIT_t), 9) = EIT_A;
Output(1:length(ECG_t), 10) = ECG_t;
Output(1:length(ECG_t), 11) = ECG_trace;
Output(1:length(PPG_t), 12) = PPG_t;
Output(1:length(PPG_t), 13) = PPG_trace0;
Output(1:length(PPG_t), 14) = PPG_trace1;
Output(1:length(PPG_t), 15) = PPG_trace2;
Output(1:length(HRB_t), 16) = HRB_t;
Output(1:length(HRB_t), 17) = HRB_HR;
Output(1:length(HRB_t), 18) = HRB_RR;
OutputTable = array2table(Output);
OutputTable.Properties.VariableNames(1:18) = {'PSD Time [s]', ...
'PSD Gauge Pressure [cmH2O]', ...
'PSD Inspiratory Differential Pressure [cmH2O]', ...
'PSD Expiratory Differential Pressure [cmH2O]', ...
'PSD Flow [L/s]', ...
'PSD Tidal Volume [L]', ...
'PSD Inspiratory Indicies' ...
...
'EIT Time [s]', ...
'EIT Global Aeration', ...
...
'ECG Time [s]', ...
'ECG [mV]', ...
...
'PPG Time [s]', ...
'PPG0', ...
'PPG1', ...
'PPG2', ...
...
'HRB Time [s]', ...
'HRB Heartrate [bpm]', ...
'HRB RR Interval [ms]'};
% writetable(OutputTable, csv_outfile_name)