Optical Pulse Oximetry [ex90x.0]¶

Calculates spO2 from given pulse oximetry signals.

../../_images/sphx_glr_fig-solution-exercise-spo2_001.png

Out:

------
Identified QRS peaks at:  [  16  130  194  256  375  471  555  636  722  809  893 1003 1117 1224
1385 1477 1584 1687 1771 1856 1947 2043 2131 2226 2319 2411 2491
2638 2698 2806 2896 2984 3076 3172 3255 3349 3462 3569 3662 3755
3935 4019 4122 4241 4330 4416 4499 4584 4665 4742 4819 4915 4997]
------
Composed signal model:  CostSegment : label: Moving Average
  └- ['Alssm : polynomial, A: (1, 1), C: (1, 1), label: None'],
  └- [Segment : a: -5, b: 5, direction: fw, g: 100, delta: 0, label: None ]
[ 127  253  373  467  552  635  719  807  892 1000 1115 1220 1306 1384
1583 1684 1769 1854 1944 2041 2130 2224 2317 2409 2489 2574 2695
2893 2982 3074 3170 3253 3346 3460 3565 3660 3754 3847 3935 4017
4238 4328 4412 4496 4583 4663 4739 4820]
[ 171  270  416  510  581  677  742  825  909 1029 1143 1238 1323 1410
1612 1702 1786 1880 1965 2066 2147 2251 2334 2425 2506 2598 2721
2911 3000 3091 3187 3270 3371 3489 3585 3677 3780 3868 3953 4040
4266 4345 4430 4523 4600 4678 4756 4838]
[80944.80545613 81006.58885978 81138.81684311 81083.55221067
83123106 80928.93968131 80824.0234319  80809.02953244
79721836 81079.05601187 81114.58230629 81147.06227972
69439485 81098.56229236 81100.19793703 81148.64704094
2839495  81015.63789444 81023.88952362 80931.47821984
47238815 80876.73942162 80903.04030187 80782.8497368
40213618 80828.6781773  80866.32157489 80985.01216733
49517527 80876.40744906 80889.20724623 80857.71657523
18821897 80924.81994765 81013.13440049 81172.50446724
59240167 81011.70031661 81089.59009723 81065.84460057
63286458 81010.63844534 81042.04401563 81044.72386864
85971366 80884.00552157 80810.94343572 80755.6640427
71868289 80944.92915366 81072.71064828]
[80521.2609977  80736.63267529 80835.70740594 80822.65397447
30562676 80630.83379386 80527.88412616 80560.15205348
46727818 80692.40793389 80788.15168944 80880.98880613
45884194 80836.85433381 80786.83443122 80779.86038152
82132016 80744.58020592 80680.65976784 80609.57145499
35632877 80594.86537942 80522.0361263  80469.63027787
45984626 80550.09733125 80524.80624123 80584.0703538
48762271 80575.90709882 80554.87865322 80558.57546765
80841118 80652.36534737 80670.95799501 80740.13455021
66479692 80714.29875048 80735.31644703 80729.8753987
01216036 80666.12974743 80582.78969642 80631.90287914
15098276 80535.77500158 80470.28838032 80435.12474312
44959734 80692.17658052 80799.22467692]
[89207.44116588 89557.12041251 89945.86336018 89872.37176049
14106287 89567.66348924 89364.38099446 89345.7140617
82477644 89992.92460197 90124.96914458 90247.12398566
89273396 90180.7705564  90195.26535718 90336.2815015
43029826 90041.20827995 90071.24328655 89881.59732423
64186834 89799.85605914 89892.98074834 89653.24194489
33592405 89725.7868757  89808.09010955 90128.56149808
08769905 89842.64332154 89874.17784167 89803.86286883
67694687 89878.53468942 90082.37759674 90450.22483724
1218083  90084.96895762 90227.31985153 90139.98740398
81383828 89971.38066195 90073.47977372 90146.06400267
76409753 89858.57034943 89711.51464235 89611.33046525
92619057 89993.39896784 90268.14967407]
[88260.91208929 88835.93483373 89181.55884028 89239.53531653
68581005 88822.67760458 88576.14111528 88654.39869064
99980461 88950.06492049 89252.94396203 89538.07389971
48103285 89484.42591921 89355.73058279 89365.03667098
02310284 89316.79061805 89192.84012951 89037.10378258
38454961 89054.09196962 88924.5705099  88821.41483646
13773695 88991.70707385 88923.44169761 89063.73170043
70226619 89047.15098189 88993.78167994 88990.73782396
52942664 89145.03879275 89170.07491432 89332.25199573
0017237  89279.38867892 89289.56063851 89272.27158415
44025176 89082.63704879 88918.88842668 89098.67351659
15519271 88954.76496104 88853.44073566 88767.29465764
63169255 89329.78832279 89537.05151854]

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
import csv as csv

import lmlib as lm
from lmlib.utils.generator import gen_slope, gen_wgn
from lmlib.utils.beta import find_max_mask
from lmlib.utils.beta import load_source_csv


y_raw = load_source_csv('./csv-data/hsp-optical_2020-11-16_08-32-48.csv') # for script start
# y_raw = load_source_csv('./csv-data/cerebral-vasoreg-diabetes-heaad-up-tilt_day1-s0030DA-noheader.csv', time_format = "H:M:S") # for script start

MIN_PEAK_DISTANCE_S = 0.6


# Signal
fs = 100
k0 = 0
K = 5000
k = range(K)
y_n = 0*gen_wgn(K, sigma=.5, seed=431412)
y_ecg = y_raw[k0:k0+K, 1]
y = (y_raw[k0:k0+K, 1:3])#  + y_n  # light intensity value


# Reference from QRS
(peaks_ecg_ind, _) = find_peaks(y_ecg, distance=fs*MIN_PEAK_DISTANCE_S)
print("------")
print("Identified QRS peaks at: ", peaks_ecg_ind)


# ALSSM models
alssm_left = lm.AlssmPoly(poly_degree=0)   # A_L, c_L
segment_left = lm.Segment(a=-5, b=5, direction=lm.FORWARD, g=100)
F = [[1]] # mixing matrix, turning on and off models per segment (1=on, 0=off)
cost = lm.CompositeCost((alssm_left, ), (segment_left,) , F, label="Moving Average")

print("------")
print("Composed signal model: ", cost)




# -- MODEL FITTING --
# Filter
separam = lm.SEParam(cost)
y = np.reshape(y, (y.shape[0], 1, -1))
separam.filter(y)

# constraint minimization
x_hat = separam.minimize_x()
y_hat = cost.eval_alssm(x_hat, alssm_selection=[1])[:,0]

y_hat_R = y_hat[:,0]
y_hat_IR = y_hat[:,1]


# find edges with conditions: only consider samples where slope increases;
(peaks_max, _) = find_peaks(y_hat_R, distance=MIN_PEAK_DISTANCE_S*fs)
peaks_max = peaks_max[1:-1] # remove first and lost to avoid border artifacts
peaks_min = np.empty_like(peaks_max)

for i, peak_max_ind in enumerate(peaks_max):
    peaks_min[i] = peak_max_ind+np.argmin(y_hat_R[peak_max_ind:int(peak_max_ind+MIN_PEAK_DISTANCE_S*fs)])



# === EXERCISE TO BE COMPLETED ==================
# SPO2 Computation
# webpage : https://omlc.org/spectra/hemoglobin/index.html
# Table : https://omlc.org/spectra/hemoglobin/summary.html


# Observed Light Intensities
# --------------------------
# at 660 nm (red)
epsilonR_Hb = 3.3  #
epsilonR_HbO2 = 0.31 #

# at 880 nm (infrared)
epsilonIR_Hb = 0.7  #
epsilonIR_HbO2 = 1.02 #  estimated


# Calibration values for R to SpO2
# --> Quadratic approximation for reference Curve MAX30101
# Source: https://www.mouser.jp/pdfDocs/SpO2-Measurement-Maxim-MAX32664-Sensor-Hub.pdf
# --------------------------
a = 1.595
b = -34.659
c = 112.689

# ======================================


R = np.log(y_hat_R[peaks_max]/y_hat_R[peaks_min]) / np.log(y_hat_IR[peaks_max]/y_hat_IR[peaks_min])

# theoretical
SaO2 = (epsilonR_Hb-epsilonIR_Hb*R) / (epsilonR_Hb- epsilonR_HbO2+ (epsilonIR_HbO2-epsilonIR_Hb)*R )

# calibrated
SaO2_Cal = 0.01 * (a*R**2 + b*R + c)


print(peaks_max)
print(peaks_min)

print(y_hat_R[peaks_max])
print(y_hat_R[peaks_min])
print(y_hat_IR[peaks_max])
print(y_hat_IR[peaks_min])


# Trajectories
trajs_max = lm.map_trajectories(cost.trajectories(x_hat[peaks_max], F=F), peaks_max, K, merge_ks=True)
trajs_min = lm.map_trajectories(cost.trajectories(x_hat[peaks_min], F=F), peaks_min, K, merge_ks=True)

# peak estimate
y_hat = cost.eval_alssm(x_hat, alssm_selection=[1])

# Plot
_, axs = plt.subplots(3, 1, sharex='all', figsize=(10, 4))

axs[0].plot(k, np.squeeze(y)[:,0], lw=1, c='black', label='')
# axs[0].plot(k, np.squeeze(y_hat)[:,0], lw=1, c='tab:gray', label='')
axs[0].plot(k, np.squeeze(trajs_max)[:,0], lw=1, c='tab:red', label='')
axs[0].plot(k, np.squeeze(trajs_min)[:,0], lw=1, c='tab:green', label='')
axs[0].scatter(peaks_max, y_hat[peaks_max,0,0], s=10, marker=7,  c='k', label='$I_{R,H}$', zorder=20)
axs[0].scatter(peaks_min, y_hat[peaks_min,0,0], s=10, marker=6,  c='k', label='$I_{R,L}$', zorder=20)

axs[0].set(xlabel='k', ylabel=r'$I-{R}$')
axs[0].set_title('Pulse Oximetry ')
axs[0].legend(loc=4)
axs[0].grid(True)


axs[1].plot(k, np.squeeze(y)[:,1], lw=1, c='black', label='$I_{IR}$')
axs[1].plot(k, np.squeeze(y_hat)[:,1], lw=1, c='tab:gray', label='')
axs[1].plot(k, np.squeeze(trajs_max)[:,1], lw=1, c='tab:red', label='')
axs[1].plot(k, np.squeeze(trajs_min)[:,1], lw=1, c='tab:green', label='')
axs[1].scatter(peaks_max, y_hat[peaks_max,0,1], s=10, marker=7,  c='k', label='$I_{IR,H}$', zorder=20)
axs[1].scatter(peaks_min, y_hat[peaks_min,0,1], s=10, marker=6,  c='k', label='$I_{R,L}$', zorder=20)

axs[1].set(xlabel='k', ylabel=r'$I-{IR}$')
axs[1].legend(loc=4)
axs[1].grid(True)


#axs[1].plot(peaks_max, R, lw=0.5, c='black', label='')
axs[2].plot(peaks_max, SaO2*100, '.--', lw=1.0, c='black', label='SpO2 (phys.)')
axs[2].plot(peaks_max, SaO2_Cal*100, '.-', lw=1.0, c='black', label='SpO2 (cal.)')
axs[2].set_ylim(80, 100)
axs[2].grid(True)
axs[2].legend(loc=4)

axs[2].set(xlabel='k', ylabel=r'$SpO2$ [%]')


plt.show()

Total running time of the script: ( 0 minutes 0.828 seconds)

Gallery generated by Sphinx-Gallery

Detecting known shapes of unknown amplitude by using polynomial approximations (Exercise 6.2) Installers