Simulating light propagation
Table of Contents
Copyright (C) 2022 Miodrag Bolic
This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details <https://www.gnu.org/licenses/>.
This code was developed by Miodrag Bolic for the book PERVASIVE CARDIAC AND RESPIRATORY MONITORING DEVICES: https://github.com/Health-Devices/CARDIAC-RESPIRATORY-MONITORING
Acknoledgement: the author would like to thank the developers of MCMatlab software.
% Changing the path from main_folder to a particular chapter
main_path=fileparts(which('Main_Content.mlx'));
if ~isempty(main_path)
%addpath(append(main_path,'/Chapter2'))
cd (append(main_path,'/Chapter6/SimulatingLightPropagation'))
addpath(append(main_path,'/Service'))
end
SAVE_FLAG=1; % saving the figures in a file
Simulating data collection for one PPG pulse
% Fig 7.10 Simulation
% The air gap of 0.01 mm
clear_all_but('SAVE_FLAG')
global zsurf;
zsurf=0.01;
Example4_BloodVessel_Pulse
for i=1:38;
model = runMonteCarlo(model);
%plotMCmatlab(model);
if model.MC.LC.res > 1
detFraction = 100*mean(mean(sum(model.MC.LC.image,3)))*model.MC.LC.fieldSize^2;
else
detFraction = 100*sum(model.MC.LC.image,3);
end
S(i)=detFraction;
end
save('SPulse_660_001_1.mat','S')
% Increasing the air gap to 0.2 mm
clear all
global zsurf;
zsurf=0.2;
Example4_BloodVessel_Pulse
for i=1:38;
model = runMonteCarlo(model);
%plotMCmatlab(model);
if model.MC.LC.res > 1
detFraction = 100*mean(mean(sum(model.MC.LC.image,3)))*model.MC.LC.fieldSize^2;
else
detFraction = 100*sum(model.MC.LC.image,3);
end
S(i)=detFraction;
end
save('SPulse_660_02_1.mat','S')
T=1/38;
t=T:T:length(S)*T;
figure
plot(t,S)
ylabel('Ilumination fraction (%)', 'FontSize', 10)
legend("660nm");
xlabel('Time (sec)', 'FontSize', 10)
% Fig 7.10 Plotting
%load("SPulse_660_02_1.mat")
load("SPulse_660_001_1.mat")
T=1/38;
t=T:T:length(S)*T;
[xData, yData] = prepareCurveData( t, S );
% Set up fittype and options.
ft = fittype( 'poly4' );
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft );
figure
plot(t,S,'o')
hold on
plot(t, fitresult(t) );
legend('MCMatlab simulation at discrete points','4th order polynomial fit')
ylabel('Ilumination fraction (%)', 'FontSize', 10)
xlabel('Time (s)', 'FontSize', 10)
title('Simulated and fitted pulse during a cardiac cycle')
annonation_save('',"Fig6.10.jpg", SAVE_FLAG);
disp('Percentage of the AC component is:')
100*(max(fitresult(t))-min(fitresult(t)))/max(fitresult(t))
ans = 1.3402
Simulation of the airgap at different frequencies
% Fig 7.7 simulation
clear all
global zsurf;
global wave_l;
global B_Cutaneous;
B=[0.2037, 0.2454];
wavelen=[660, 940];
depth=[0.001, 0.01, 0.2, 0.3]; % 0.4 is too much
for j=1:2
for k=1:2
for i=1:4;
zsurf=depth(i);
wave_l=wavelen(j);
B_Cutaneous=B(k);
clear_all_but('zsurf','wave_l','B_Cutaneous', 'i', 'j', 'k', 'S','depth','wavelen','B');
Example4_BloodVessel_changeS_barrier1
model = runMonteCarlo(model);
plotMCmatlab(model);
if model.MC.LC.res > 1
detFraction = 100*mean(mean(sum(model.MC.LC.image,3)))*model.MC.LC.fieldSize^2;
else
detFraction = 100*sum(model.MC.LC.image,3);
end
S(j,k,i)=detFraction;
end
end
end
save('SBarrier1.mat','S')
% Fig 7.7 plotting
load('SBarrier1.mat')
depth=[0.001, 0.01, 0.2, 0.3]; % 0.4 is too much
figure
N=3;
Dia660=reshape(S(1,1,:),[4,1]);
Sys660=reshape(S(1,2,:),[4,1]);
Dia940=reshape(S(2,1,:),[4,1]);
Sys940=reshape(S(2,2,:),[4,1]);
plot(depth(1:N),Dia660(1:N))
hold on
plot(depth(1:N),Sys660(1:N))
plot(depth(1:N),Dia940(1:N))
plot(depth(1:N),Sys940(1:N))
legend('660 nm diastolic point','660 nm systolic point','940 nm diastolic point','940 nm systolic point')
ylabel('Ilumination fraction (%)', 'FontSize', 10)
xlabel('Barrier depth (mm)', 'FontSize', 10)
figure
plot(depth(1:N),((Dia660(1:N)-Sys660(1:N))./Dia660(1:N))./((Dia940(1:N)-Sys940(1:N))./Dia940(1:N)),'o-')
ylabel('Ratio of ratios', 'FontSize', 10)
xlabel('Air gap depth (mm)', 'FontSize', 10)
title('Ratio of ratios vs. the air gap depth')
annonation_save('b)',"Fig6.7b.jpg", SAVE_FLAG);
figure
plot(depth(1:N),100*((Dia660(1:N)-Sys660(1:N))./Dia660(1:N)),'o-')
hold on
plot(depth(1:N),100*((Dia940(1:N)-Sys940(1:N))./Dia940(1:N)),'*-')
legend('660 nm','940 nm')
ylabel('Ratio AC over DC (%)', 'FontSize', 10)
xlabel('Air gap depth (mm)', 'FontSize', 10)
title('AC over DC ratios for red and infrared light')
annonation_save('a)',"Fig6.7a.jpg", SAVE_FLAG);
Plotting the normalized absorbed power and normalized fluence rate of collected light
%Fig 7.3 and 7.8
clear_all_but('SAVE_FLAG')
global zsurf;
zsurf=0.01;
Example4_BloodVessel_Pulse
model = runMonteCarlo(model);
plotMCmatlab(model);
if model.MC.LC.res > 1
detFraction = 100*mean(mean(sum(model.MC.LC.image,3)))*model.MC.LC.fieldSize^2;
else
detFraction = 100*sum(model.MC.LC.image,3);
end
S(i)=detFraction;
Effects of air gap to the ratio of ratios
% Fig 7.9 simulation
clear all
global ii;
global wave_l;
global B_Cutaneous;
B=[0.2037, 0.2454];
wavelen=[660, 940];
for j=1:2
for k=1:2
clear_all_but('wave_l','B_Cutaneous', 'i1', 'j', 'k', 'S','B','wavelen','ii');
wave_l=wavelen(j);
B_Cutaneous=B(k);
ii=0.45;
Example4_BloodVessel_changeS_new1
for i1=1:11;
ii=0.45+0.05*i1;
model = runMonteCarlo(model);
%plotMCmatlab(model);
if model.MC.LC.res > 1
detFraction = 100*mean(mean(sum(model.MC.LC.image,3)))*model.MC.LC.fieldSize^2;
else
detFraction = 100*sum(model.MC.LC.image,3);
end
S(j,k,i1)=detFraction;
end
end
end
save('SChange.mat','S')
test=1;
% Fig 7.9 plotting
load('SChange.mat')
SpO2=50:5:100;
RatioOfRatios=((reshape(S(1,2,:),[11,1])-reshape(S(1,1,:),[11,1]))./reshape(S(1,2,:),[11,1]))./((reshape(S(2,2,:),[11,1])-reshape(S(2,1,:),[11,1]))./reshape(S(2,2,:),[11,1]));
%% Fit: 'untitled fit 1'.
p = polyfit(RatioOfRatios',SpO2,1);
y1 = polyval(p,RatioOfRatios');
figure
plot(RatioOfRatios,SpO2,'o')
hold on
plot(RatioOfRatios,y1,'-')
xlabel('Ratio of ratios', 'FontSize', 10)
ylabel('Saturation (%)', 'FontSize', 10)
legend('SpO_2 estimate','SpO_2=C_B*R+C_A')
ylim([50,100])
title('Estimated and fitted SpO_2 vs. the ratio of ratios')
annonation_save('b)',"Fig6.9b.jpg", SAVE_FLAG);
figure
plot(SpO2,reshape(S(1,1,:),[11,1]),'o-')
hold on
plot(SpO2,reshape(S(1,2,:),[11,1]),'*-')
plot(SpO2,reshape(S(2,1,:),[11,1]),'o-')
plot(SpO2,reshape(S(2,2,:),[11,1]),'*-')
legend('660 nm diastolic','660 nm systolic','940 nm diastolic','940 nm systolic',"Location","best")
ylabel('Ilumination fraction (%)', 'FontSize', 10)
xlabel('Saturation (%)', 'FontSize', 10)
title('Illumination vs SpO_2 for systolic and diastolic points')
annonation_save('a)',"Fig6.9a.jpg", SAVE_FLAG);
