In [1]:
import math
%matplotlib inline
import matplotlib
import matplotlib.pyplot as pp
import numpy as np
import scipy.signal
In [2]:
Fs = 2000
Fs_rsp = 90
Fs_ppg = 200

F_rsp_target = 1/15

file = "ppg-rsp-data/ppg-rsp-resonate-15-180sec.txt"
data = np.genfromtxt(fname=file, delimiter=",").T
ppg = scipy.signal.resample_poly(data[0], Fs_ppg, Fs)
rsp = scipy.signal.resample_poly(data[1], Fs_rsp, Fs)
In [3]:
pp.plot(np.linspace(0, 1, len(rsp)), rsp)
pp.plot(np.linspace(0, 1, len(ppg)), 40*ppg)
Out[3]:
[<matplotlib.lines.Line2D at 0x7f788be8aba8>]
In [208]:
b, a = scipy.signal.butter(2, [1/60, 1/10], btype = "bandpass", fs = Fs_rsp)
rsp_filtered = scipy.signal.lfilter(b, a, rsp)

pp.axhline(0)
pp.plot(rsp, color = "lightgrey")
pp.plot(rsp_filtered)

(b, a)
Out[208]:
(array([ 6.82859420e-06,  0.00000000e+00, -1.36571884e-05,  0.00000000e+00,
         6.82859420e-06]),
 array([ 1.        , -3.99258209,  5.97778678, -3.97782723,  0.99262254]))
In [224]:
# TODO this is based on the frequency
gain = Fs_rsp/(2*math.pi*F_rsp_target)

deriv = scipy.signal.lfilter([gain/2, 0, -gain/2], [1], rsp_filtered)
# deriv = scipy.signal.lfilter([gain, -gain], [1], rsp_filtered)

pp.figure(figsize = (20, 10))
pp.axhline(0)
pp.plot(rsp_filtered)
pp.plot(deriv)

RSP_MAG = np.sqrt(rsp_filtered*rsp_filtered + deriv*deriv)
b, a = scipy.signal.butter(2, 1/60, btype = "lowpass", fs = Fs_rsp)
RSP_MAG_FILTERED = scipy.signal.lfilter(b, a, RSP_MAG)

pp.plot(RSP_MAG, color = "pink")
pp.plot(RSP_MAG_FILTERED, color = "red")

(a, b)
Out[224]:
(array([ 1.        , -1.99851904,  0.99852014]),
 array([2.73952796e-07, 5.47905592e-07, 2.73952796e-07]))
In [225]:
rsp_peak_x, _ = scipy.signal.find_peaks(rsp_filtered, width = Fs_rsp)
rsp_peaks = np.array([rsp_peak_x, rsp_filtered.take(rsp_peak_x)])
pp.plot(rsp_filtered)
pp.plot(rsp_peaks[0], rsp_peaks[1], "o")
(len(rsp)/len(rsp_peak_x))/Fs_rsp
Out[225]:
7.865714285714286
In [229]:
LEN = len(rsp)
dt = (2*math.pi*F_rsp_target)/Fs_rsp
# dt /= 1.1
# dt *= 1.1
c = 0.2
damping = math.exp(-c*dt)

sin, cos = 0, 0
SIN, COS = np.empty(LEN), np.empty(LEN)
SIN[0], COS[0] = sin, cos
UNWARPED = np.empty(LEN)
UNWARPED[0] = 0

for i in range(1, LEN):
    mag = math.sqrt(cos*cos + sin*sin) + 1e-5
    warp = 1 - 0.2*sin/mag
    rsp_sin = rsp_filtered[i]*warp
    UNWARPED[i] = rsp_sin
    
    cos *= damping
    cos -= dt*(sin + c*rsp_sin)
    sin *= damping
    sin += dt*(cos + c*deriv[i])
    SIN[i], COS[i] = sin, cos

pp.figure(figsize = (20,5))
pp.axhline(0, color = "lightgrey")

X = np.linspace(0, dt*LEN, LEN)
pp.plot(X, UNWARPED, color = "teal")
# pp.plot(X, deriv, color = "orange")
# pp.plot(X, SIN, color = "lightgrey")
pp.plot(X, -COS, color = "darkgrey")

SHM_MAG = np.sqrt(SIN*SIN + COS*COS) + 1e-50
# pp.plot(X, RSP_MAG_FILTERED, color = "green")
# pp.plot(X, SHM_MAG, color = "red")
pp.axhline(1, color = "green")
pp.plot(X, SHM_MAG/(RSP_MAG_FILTERED + 1e-1), color = "red")

{
    "dt": dt,
    "1/F_rsp_target": 1/F_rsp_target,
    "avg": np.average(SHM_MAG/(RSP_MAG_FILTERED + 1e-1)),
}
Out[229]:
{'dt': 0.0041887902047863905,
 '1/F_rsp_target': 15.0,
 'avg': 0.6848355742348389}
In [230]:
pp.figure(figsize = (10, 10))
pp.plot(rsp_filtered, deriv, color = "lightblue")
pp.plot(COS, SIN, color = "orange")
pp.axis("equal")

pp.plot(UNWARPED, deriv)
Out[230]:
[<matplotlib.lines.Line2D at 0x7f788f7b9c18>]
In [9]:
b, a = scipy.signal.butter(2, [0.5, 4], btype = "bandpass", fs = Fs_ppg)
ppg_filtered = scipy.signal.lfilter(b, a, ppg)

r = range(0, len(ppg))
pp.plot(ppg[r])
pp.plot(ppg_filtered[r])
Out[9]:
[<matplotlib.lines.Line2D at 0x7f78959307b8>]