%% Data Processing for April 2023 Trial
% Ella Guy
% Last Updated: 29AUG2023
clear all, close all, clc
% Load Data
%--------------------------------------------------------------------------
% Inputs-------------------------------------------------------------------
for SubjectNumber = [1:68,70:80] % de-identifyed subject number
% sensor offsets
OffsetA = -0.150;
OffsetB = -0.064;
% processing criteria values
MinProm = 0.05; % 50 mL
MeanThreshold = 5000;
Tapethreshold = 1;
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
EIT_infile_format = 'EIT_rawData/S%02d_01_001_01.bin';
EIT_infile = sprintf(EIT_infile_format, SubjectNumber);
EIT_data = read_binData(EIT_infile);
% EIT_Infile:
BOB_infile_format = 'PQ_rawData/Subject%d.mat';
BOB_infile = sprintf(BOB_infile_format, SubjectNumber);
load(BOB_infile)
% BOB_Infile:
% 'time' - time [s]
% 'GaugeP' - Gauge pressure reading [cmH2O]
% 'InhaleDeltaP' - Inspiratory differenrtial pressure [cmH2O]
% 'ExhaleDeltaP' - Expiratory differenrtial pressure [cmH2O]
% 'Chest' - Chest circumference [mm]
% 'Abd' - Abdominal circumference [mm]
% 'Depth_chest' - Chest crossectional depth (a) [cm]
% 'Width_chest' - Chest crossectional width (b) [cm]
% EIT Data Processing -----------------------------------------------------
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)]);
% BOB Data Processing -----------------------------------------------------
BOB_sample_frequency = 100; % [Hz]
BOB_timestep = 1/BOB_sample_frequency; % [s]
% Trimming
for ai = 1:length(time)
if time(ai) < 0
TE = ai-1;
break
end
end
GaugeP = GaugeP(1:TE);
InhaleDeltaP = InhaleDeltaP(1:TE) - OffsetA;
ExhaleDeltaP = ExhaleDeltaP(1:TE) - OffsetB;
Chest = Chest(1:TE);
Abd = Abd(1:TE);
Time = BOB_timestep.*([1:length(GaugeP)]) - BOB_timestep;
% 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 = ExhaleDeltaP - InhaleDeltaP;
Flow = zeros(length(DeltaP), 1);
for bi = 1 : length(DeltaP)
if DeltaP(bi) < 0
Flow(bi) = 1000*cd*A2*sqrt((2*(-1.*DeltaP(bi)*98.0665))/(rho*(1-(A2/A1)^2)));
else
Flow(bi) = -1000*cd*A2*sqrt((2*(DeltaP(bi)*98.0665))/(rho*(1-(A2/A1)^2)));
end
end
% Indicies
VC = cumtrapz(BOB_timestep, Flow);
VCM = movmean(VC,MeanThreshold);
CriterionA = islocalmin(VC-VCM, 'MinProminence', MinProm, 'FlatSelection', 'last');
InspInd = [];
for bii = 2:length(Flow)
if CriterionA(bii) == 1
InspInd = [InspInd, bii];
end
end
% % troubleshooting plots -------------------------------------------------
figure(7)
hold on
plot(VC)
plot(VCM)
plot(InspInd, VC(InspInd), '.')
ylabel("Cumulative Volume")
hold off
%
% figure(8)
% hold on
% plot(VC-VCM)
% plot(InspInd, VC(InspInd)-VCM(InspInd), '.')
% ylabel("Cumulative Volume Normalised By Moving Mean")
% hold off
% % -----------------------------------------------------------------------
% Breathwise Tidal Volume
V_working = 0;
Volume = [V_working];
for biv = 2:length(Flow)
if ismember(biv, InspInd) == 0
V_working = V_working + trapz(BOB_timestep, Flow(biv-1:biv));
else
V_working = 0;
end
Volume = [Volume; V_working];
end
% Plot axes variables
ChestMin = floor(max(Chest(1:900)));
AbdominalMax = ceil(max(Abd(1:900)));
% PEEP Profile ------------------------------------------------------------
PEEP_Profile = [0,0, 4,4, 4.5,4.5, 5,5, 5.5,5.5, 6,6, ...
6.5,6.5, 7,7, 7.5,7.5, 8,8, 8.5,8.5, 9,9, 9.5,9.5, ...
10,10, 10.5,10.5, 11,11, 11.5,11.5, 12,12, 0,0];
PEEP_Timings = [0,60, 90,120, 140,170, 190,220, 240,270, 290,310, ...
340,370, 390,420, 440,470, 490,520, 540,570, 590,620, 640,670, ...
690,720, 740,770, 790,820, 840,870, 890,920, 940,1000];
% Plot Data #1 (Pre-Alignment)---------------------------------------------
figure(1)
ax(1) = subplot(6, 1, 1);
hold on
plot(Time, GaugeP)
plot(Time(InspInd), GaugeP(InspInd), 'x')
hold off
grid on
grid minor
xlim([0 1200])
ylim([-20 20])
title("Gauge Pressure [cmH_2O]", 'Fontsize', 12)
ax(2) = subplot(6, 1, 2);
hold on
plot(Time, Flow)
hold off
grid on
grid minor
xlim([0 1200])
title("Flow [L/s]", 'Fontsize', 12)
ax(3) = subplot(6, 1, 3);
hold on
plot(Time, Volume)
plot(Time(InspInd), Volume(InspInd), 'x')
hold off
grid on
grid minor
ylim([0 3])
xlim([0 1200])
title("Volume [L]", 'Fontsize', 12)
ax(4) = subplot(6, 1, 4);
plot(Time, Chest)
grid on
grid minor
xlim([0 1200])
ylim([ChestMin-10 ChestMin+30])
title("Chest Circumference [mm]", 'Fontsize', 12)
ax(5) = subplot(6, 1, 5);
plot(Time, Abd)
grid on
grid minor
xlim([0 1200])
ylim([AbdominalMax-30 AbdominalMax+10])
title("Abdominal Circumference [mm]", 'Fontsize', 12)
ax(6) = subplot(6, 1, 6);
plot(Time_Aeration, Aeration)
grid on
grid minor
xlim([0 1200])
title("Global Aeration (EIT)", 'Fontsize', 12)
linkaxes(ax, 'x')
%--------------------------------------------------------------------------
% % troubleshooting plots -------------------------------------------------
% figure(2)
%
% hold on
% plot(Time, GaugeP)
% plot(PEEP_Timings, 0.5*(PEEP_Profile))
% ylabel("Pressure Profile comparison")
% hold off
% grid on
% grid minor
% xlim([0 1200])
% ylim([-20 20])
%--------------------------------------------------------------------------
% alignment working plots -------------------------------------------------
figure(3)
ax1(1) = subplot(2, 1, 1);
hold on
plot(Time, Volume)
plot(Time(InspInd), Volume(InspInd), 'x')
title("Flow-based Volume")
hold off
grid on
grid minor
xlim([0 1200])
ylim([-2 2])
ax1(2) = subplot(2, 1, 2);
plot(Time_Aeration, Aeration)
title("Aeration")
grid on
grid minor
xlim([0 1200])
linkaxes(ax1, 'x')
%--------------------------------------------------------------------------
% Dynamic Circumference Processing
countC = 0;
countA = 0;
for a = 2:length(Chest)
chestjump = Chest(a)- Chest(a-1);
if abs(chestjump) > Tapethreshold
Chest(a:end) = Chest(a:end) - chestjump;
countC = countC + 1;
end
abdjump = Abd(a)-Abd(a-1);
if abs(abdjump) > Tapethreshold
Abd(a:end) = Abd(a:end) - abdjump;
countA = countA + 1;
end
end
countC
countA
%% Aligning Recordings
% Align -------------------------------------------------------------------
S01 = [188,189,188,189];
S02 = [631,632,630,632];
S03 = [889,890,887,888];
S04 = [50,51,50,51];
S05 = [293,294,293,294];
S06 = [178,179,177,179];
S07 = [6,7,2,3];
S08 = [189,190,188,189];
S09 = [335,336,335,336];
S10 = [731,732,730,731];
S11 = [90,91,89,90];
S12 = [514,515,513,514];
S13 = [44,46,52,53];
S14 = [757,758,756,757];
S15 = [662,663,662,663];
S16 = [693,696,692,695];
S17 = [106,107,106,107];
S18 = [293,294,292,293];
S19 = [458,459,475,476];
S20 = [81,82,80,81];
S21 = [53,54,51,52];
S22 = [66,67,65,66];
S23 = [923,924,923,924];
S24 = [529,530,528,529];
S25 = [329,330,328,329];
S26 = [337,339,337,339];
S27 = [75,76,76,77];
S28 = [23,25,23,24];
S29 = [296,297,296,297];
S30 = [333,334,332,333];
S31 = [66,68,65,66];
S32 = [81,82,81,82];
S33 = [256,257,255,256];
S34 = [195,196,208,209];
S35 = [877,878,875,876];
S36 = [65,66,64,65];
S37 = [381,382,381,382];
S38 = [243,246,243,244];
S39 = [448,450,447,449];
S40 = [507,508,506,507];
S41 = [45,46,44,45];
S42 = [573,573.8,571,573];
S43 = [640,641,639,640];
S44 = [184,186,183,185];
S45 = [829,830,827,828];
S46 = [477,478,476,477];
S47 = [254,256,255,257];
S48 = [143,147,143,145];
S49 = [466,468,465,467];
S50 = [74,75,72,74];
S51 = [759,780,758,759];
S52 = [515,516,514,515];
S53 = [690,692,691,693];
S54 = [322,323,322,323];
S55 = [724,725,723,724];
S56 = [596,599,596,598];
S57 = [78,79,77,78];
S58 = [71,73,71,73];
S59 = [539,541,539,541];
S60 = [720,721,718,719];
S61 = [6.5,7.5,6.5,7.5];
S62 = [1.5,4,0.5,2];
S63 = [977,979,976,978];
S64 = [58,60,58,60];
S65 = [12,15,12,15];
S66 = [761,762,759,760];
S67 = [958,960,978,980];
S68 = [114,116,113,115];
S69 = [0,0,0,0];
S70 = [67,68,66,67];
S71 = [960,961,966,967];
S72 = [916,917,929,930];
S73 = [750,751,749,750];
S74 = [58,61,53,56];
S75 = [45,50,45,50];
S76 = [1026,1030,1026,1030];
S77 = [207,209,207,209];
S78 = [160,162,160,162];
S79 = [187,188,188,189];
S80 = [924,925,924,925];
Align = [S01; S02; S03; S04; S05; S06; S07; S08; S09; S10; ...
S11; S12; S13; S14; S15; S16; S17; S18; S19; S20; ...
S21; S22; S23; S24; S25; S26; S27; S28; S29; S30; ...
S31; S32; S33; S34; S35; S36; S37; S38; S39; S40; ...
S41; S42; S43; S44; S45; S46; S47; S48; S49; S50; ...
S51; S52; S53; S54; S55; S56; S57; S58; S59; S60; ...
S61; S62; S63; S64; S65; S66; S67; S68; S69; S70; ...
S71; S72; S73; S74; S75; S76; S77; S78; S79; S80];
%--------------------------------------------------------------------------
% Alignment Breath --------------------------------------------------------
% Range input for last Peak Vol before unplugging CPAP
% OR// 1st Peak Vol at 4cmh2O
% OR// other distinctive breath
% In BOB data
T_min_BOB = Align(SubjectNumber, 1); %[s]
T_max_BOB = Align(SubjectNumber, 2); %[s]
% In EIT data
T_min_EIT = Align(SubjectNumber, 3); %[s]
T_max_EIT = Align(SubjectNumber, 4); %[s]
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Range of first peak at 4cmH2O PEEP
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(Volume(I_min_BOB:I_max_BOB));
[LM_EIT, LMI_EIT] = max(Aeration(I_min_EIT:I_max_EIT));
LMT_BOB = Time(I_min_BOB + LMI_BOB);
LMT_EIT = Time_Aeration(I_min_EIT + LMI_EIT);
% Aligning
T_align = LMT_EIT - LMT_BOB;
if T_align < 0
Time = Time' + T_align;
I_align = round(abs(T_align*BOB_sample_frequency));
T = Time(I_align:end);
P = GaugeP(I_align:end);
Q = Flow(I_align:end);
V = Volume(I_align:end);
CC = Chest(I_align:end);
AC = Abd(I_align:end);
Ind_working = InspInd' - I_align;
for ci = 1:length(Ind_working)
if Ind_working(ci) >= 0
break
end
end
Ind = Ind_working(ci:end);
T_GA = Time_Aeration;
GA = Aeration;
Ind_GA_working = floor(Ind./2);
for cii = 1:length(Ind_GA_working)
if Ind_GA_working(cii) > length(GA)
break
end
end
Ind_GA = Ind_GA_working(1:cii-1);
else
I_align = round(abs(T_align*EIT_sample_frequency)+1);
T = Time;
P = GaugeP;
Q = Flow;
V = Volume;
CC = Chest;
AC = Abd;
Ind_working = InspInd';
T_GA = Time_Aeration(I_align:end) - Time_Aeration(I_align);
GA = Aeration(I_align:end);
Ind_GA_working = floor(Ind./2);
for cii = 1:length(Ind_GA_working)
if Ind_GA_working(cii) > length(GA)
break
end
end
Ind_GA = Ind_GA_working(1:cii-1);
end
for ci = 1:length(Ind_working)
if Ind_working(ci) > 0
break
end
end
Ind = Ind_working(ci:end);
Ind_GA = Ind_GA(ci:end);
% Plot Data #2 (Alignment Check)-------------------------------------------
% ***Check PEEP profile against Gauge Pressure for cueing alignment offset
figure(4)
ax(1) = subplot(6, 1, 1);
hold on
plot(T, P)
plot(T(Ind), P(Ind), 'x')
hold off
grid on
grid minor
xlim([0 1200])
ylim([-20 20])
title("Gauge Pressure [cmH_2O]", 'Fontsize', 12)
ax(2) = subplot(6, 1, 2);
hold on
plot(T, Q)
plot(T(Ind), Q(Ind), 'x')
hold off
grid on
grid minor
xlim([0 1200])
title("Flow [L/s]", 'Fontsize', 12)
ax(3) = subplot(6, 1, 3);
hold on
plot(T, V)
plot(T(Ind), V(Ind), 'x')
hold off
grid on
grid minor
xlim([0 1200])
title("Volume [L]", 'Fontsize', 12)
ax(4) = subplot(6, 1, 4);
hold on
plot(T, CC)
plot(T(Ind), CC(Ind), 'x')
hold off
grid on
grid minor
xlim([0 1200])
ylim([ChestMin-10 ChestMin+30])
title("Chest Circumference [mm]", 'Fontsize', 12)
ax(5) = subplot(6, 1, 5);
hold on
plot(T, AC)
plot(T(Ind), AC(Ind), 'x')
hold off
grid on
grid minor
xlim([0 1200])
ylim([AbdominalMax-30 AbdominalMax+10])
title("Abdominal Circumference [mm]", 'Fontsize', 12)
ax(6) = subplot(6, 1, 6);
hold on
plot(T_GA, GA)
plot(T_GA(Ind_GA), GA(Ind_GA), 'x')
hold off
grid on
grid minor
xlim([0 1200])
title("Global Aeration (EIT)", 'Fontsize', 12)
linkaxes(ax, 'x')
% alignment working plots -------------------------------------------------
figure(5)
hold on
plot(T, P)
plot(PEEP_Timings, 0.5*(PEEP_Profile))
title("Pressure Profile Comparison")
hold off
grid on
grid minor
xlim([0 1200])
ylim([-20 20])
%--------------------------------------------------------------------------
% % troubleshooting plots -------------------------------------------------
figure(5)
ax(1) = subplot(3, 1, 1);
hold on
plot(Time, InhaleDeltaP)
title("Inhale Differential Pressure")
hold off
grid on
grid minor
ax(2) = subplot(3, 1, 2);
hold on
plot(Time, ExhaleDeltaP)
title("Exhale Differential Pressure")
hold off
grid on
grid minor
ax(2) = subplot(3, 1, 3);
hold on
plot(Time, DeltaP)
title("Combined Differential Pressure")
hold off
grid on
grid minor
%% Saves Data
outfile_format = 'ProcessedDataset/ProcessedData_Subject%02d_Indexed.mat';
outfile_name = sprintf(outfile_format, SubjectNumber);
save(outfile_name, 'T', 'P', 'Q', 'V', 'CC', 'AC', 'Ind', ...
'T_GA', 'GA', 'Ind_GA')
csv_outfile_format = 'ProcessedDataset/ProcessedData_Subject%02d.csv';
csv_outfile_name = sprintf(csv_outfile_format, SubjectNumber);
Output = zeros(length(T), 10);
Output(:,1) = T;
Output(:,2) = P;
Output(:,3) = Q;
Output(:,4) = V;
Output(:,5) = CC;
Output(:,6) = AC;
Output(1:length(Ind), 7) = Ind;
Output(1:length(T_GA), 8) = T_GA;
Output(1:length(GA), 9) = GA;
Output(1:length(Ind_GA), 10) = Ind_GA;
OutputTable = array2table(Output);
OutputTable.Properties.VariableNames(1:10) = {'Time [s]', ...
'Pressure [cmH2O]', 'Flow [L/s]', 'V_tidal [L]', ...
'Chest [mm]', 'Abd [mm]', 'Inspiratory Indicies' ...
'Time (Aeration Data)_[s]', 'Global Aeration', ...
'Inspiratory Indicies (Aeration Data)'};
writetable(OutputTable, csv_outfile_name)
writetable(OutputTable, csv_outfile_name)
end