%% 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