ABP notch localization [ex905.0]

Identifies the notch in the arterial blood pressure (ABP) signal

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

Out:

------
Identified QRS peaks at:  [  16  130  194  256  375  471  555  636  722  809  893 1003 1117 1224
 1309 1385 1477 1584 1687 1771 1856 1947 2043 2131 2226 2319 2411 2491
 2578 2638 2698 2806 2896 2984 3076 3172 3255 3349 3462 3569 3662 3755
 3850 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
 1475 1583 1684 1769 1854 1944 2041 2130 2224 2317 2409 2489 2574 2695
 2802 2893 2982 3074 3170 3253 3346 3460 3565 3660 3754 3847 3935 4017
 4120 4238 4328 4412 4496 4583 4663 4739 4820]
[ 171  270  416  510  581  677  742  825  909 1029 1143 1238 1323 1410
 1499 1612 1702 1786 1880 1965 2066 2147 2251 2334 2425 2506 2598 2721
 2821 2911 3000 3091 3187 3270 3371 3489 3585 3677 3780 3868 3953 4040
 4148 4266 4345 4430 4523 4600 4678 4756 4838]
[80944.80545613 81006.58885978 81138.81684311 81083.55221067
 81058.83123106 80928.93968131 80824.0234319  80809.02953244
 80888.79721836 81079.05601187 81114.58230629 81147.06227972
 81117.69439485 81098.56229236 81100.19793703 81148.64704094
 81084.2839495  81015.63789444 81023.88952362 80931.47821984
 80959.47238815 80876.73942162 80903.04030187 80782.8497368
 80799.40213618 80828.6781773  80866.32157489 80985.01216733
 80941.49517527 80876.40744906 80889.20724623 80857.71657523
 80976.18821897 80924.81994765 81013.13440049 81172.50446724
 81062.59240167 81011.70031661 81089.59009723 81065.84460057
 81023.63286458 81010.63844534 81042.04401563 81044.72386864
 80901.85971366 80884.00552157 80810.94343572 80755.6640427
 80798.71868289 80944.92915366 81072.71064828]
[80521.2609977  80736.63267529 80835.70740594 80822.65397447
 80754.30562676 80630.83379386 80527.88412616 80560.15205348
 80631.46727818 80692.40793389 80788.15168944 80880.98880613
 80899.45884194 80836.85433381 80786.83443122 80779.86038152
 80787.82132016 80744.58020592 80680.65976784 80609.57145499
 80622.35632877 80594.86537942 80522.0361263  80469.63027787
 80516.45984626 80550.09733125 80524.80624123 80584.0703538
 80599.48762271 80575.90709882 80554.87865322 80558.57546765
 80667.80841118 80652.36534737 80670.95799501 80740.13455021
 80728.66479692 80714.29875048 80735.31644703 80729.8753987
 80729.01216036 80666.12974743 80582.78969642 80631.90287914
 80575.15098276 80535.77500158 80470.28838032 80435.12474312
 80574.44959734 80692.17658052 80799.22467692]
[89207.44116588 89557.12041251 89945.86336018 89872.37176049
 89850.14106287 89567.66348924 89364.38099446 89345.7140617
 89491.82477644 89992.92460197 90124.96914458 90247.12398566
 90191.89273396 90180.7705564  90195.26535718 90336.2815015
 90183.43029826 90041.20827995 90071.24328655 89881.59732423
 89974.64186834 89799.85605914 89892.98074834 89653.24194489
 89692.33592405 89725.7868757  89808.09010955 90128.56149808
 89989.08769905 89842.64332154 89874.17784167 89803.86286883
 90029.67694687 89878.53468942 90082.37759674 90450.22483724
 90187.1218083  90084.96895762 90227.31985153 90139.98740398
 90015.81383828 89971.38066195 90073.47977372 90146.06400267
 89834.76409753 89858.57034943 89711.51464235 89611.33046525
 89677.92619057 89993.39896784 90268.14967407]
[88260.91208929 88835.93483373 89181.55884028 89239.53531653
 89066.68581005 88822.67760458 88576.14111528 88654.39869064
 88799.99980461 88950.06492049 89252.94396203 89538.07389971
 89610.48103285 89484.42591921 89355.73058279 89365.03667098
 89401.02310284 89316.79061805 89192.84012951 89037.10378258
 89107.38454961 89054.09196962 88924.5705099  88821.41483646
 88946.13773695 88991.70707385 88923.44169761 89063.73170043
 89096.70226619 89047.15098189 88993.78167994 88990.73782396
 89199.52942664 89145.03879275 89170.07491432 89332.25199573
 89310.0017237  89279.38867892 89289.56063851 89272.27158415
 89243.44025176 89082.63704879 88918.88842668 89098.67351659
 88989.15519271 88954.76496104 88853.44073566 88767.29465764
 89066.63169255 89329.78832279 89537.05151854]
[0.49181478 0.41285366 0.43857742 0.45607807 0.42977858 0.44183344
 0.41432015 0.39710671 0.41057895 0.4101084  0.41473798 0.4163768
 0.41655745 0.41697787 0.41398334 0.42137712 0.42036987 0.41487739
 0.43317089 0.42218613 0.43081239 0.41865271 0.43581891 0.4167573
 0.41989386 0.42026828 0.42752167 0.41759761 0.42487256 0.41855054
 0.42072944 0.40749595 0.41188624 0.41155191 0.41582043 0.42942395
 0.42236981 0.40943814 0.41908699 0.42934289 0.42273075 0.42929639
 0.44050062 0.43696668 0.42786174 0.42680749 0.43954146 0.42026221
 0.40636754 0.42254795 0.41551631]
------

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]

# Square Error and LCR
#error_edge = separam.eval_errors(x_hat_edge)
#lcr = -1/2 * np.log(np.divide(error_edge, error_line))

# 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)])
#(peaks_min, _) = find_peaks(-y_hat_R, distance=MIN_PEAK_DISTANCE_S*fs)



# SPO2 Computation

# webpage : https://omlc.org/spectra/hemoglobin/index.html
# Tabelle hier : https://omlc.org/spectra/hemoglobin/summary.html


epsilonR_Hb = 3.3# 0.68  # 660nm (rot)
epsilonR_HbO2 = 0.31 # geschätzter Wert damit Ergebnis stimmt # 319# 0.09 # 0.073

epsilonIR_Hb = 0.7  # 0.19 #  880nm
epsilonIR_HbO2 = 1.02# 0.24 # 0.26 # 0.29

# calibration values for R to SpO2
# 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])
#R = (y_hat_R[peaks_max] - y_hat_R[peaks_min]) / (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])

print(R)

print("------")

# 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.708 seconds)

Gallery generated by Sphinx-Gallery