Fitting of polynomial model of different orders to ECG signal [ex90x.0]

../../_images/sphx_glr_fig-solution-exercise-model-baseline_001.png

Out:

CostSegment : label: None
  └- ['Alssm : polynomial, A: (3, 3), C: (1, 3), label: alssm-pulse', 'Alssm : polynomial, A: (3, 3), C: (1, 3), label: alssm-baseline'],
  └- [Segment : a: -inf, b: -1, direction: fw, g: 100, delta: 0, label: segment 1 , Segment : a: 0, b: 80, direction: fw, g: 2000, delta: 0, label: segment 2 , Segment : a: 81, b: inf, direction: bw, g: 100, delta: 81, label: segment 3 ]
H_A :  [[0. 0. 0. 0. 0.]
 [1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]
H_0 :  [[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]
 [1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

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

import lmlib as lm
from lmlib.utils.generator import *

# --------------- loading test signal -----------------------
file_name = 'EECG_BASELINE_9CH_10S_FS2400HZ.csv'
K = 10000  # number of samples to process
y = load_multi_channel(file_name, K, ch_select=[3, ])
y = y[0:]

# --------------- parameters of example -----------------------
LCR_THD = 0.01  # minimum log-cost ratio to detect a pulse in noise

SHAPE_LEN = 80  # length of shape to be found
MIN_DIST = 1500  # minimal number of samples between two pulses

g_sp = 2000  # pulse window weight, effective sample number under the window # (larger value lead to a more rectangular-like windows while too large values might lead to nummerical instabilities in the recursive computations.)
g_bl = 100  # baseline window weight, effective sample number under the window (larger value leads to a wider window)

N1 = 3  # number of polynomial coefficient for baseline  Note: N>4 might leads to nummerical instabilities in the recursions
N2 = 3  # number of polynomial coefficient for spike. Note: N>4 might leads to  nummerical instabilities in the recursions.

# --------------- main -----------------------

# Defining ALSSM models
alssm_baseline = lm.AlssmPoly(poly_degree=N1 - 1, label="alssm-baseline")
alssm_pulse = lm.AlssmPoly(poly_degree=N2 - 1, label="alssm-pulse")

# Defining segments with a left- resp. right-sided decaying window and a center segment with nearly rectangular window
segmentL = lm.Segment(a=-np.infty, b=-1, direction=lm.FORWARD, g=g_bl, delta=0,
                      label='segment 1')
segmentC = lm.Segment(a=0, b=SHAPE_LEN, direction=lm.FORWARD, g=g_sp, label='segment 2')
segmentR = lm.Segment(a=SHAPE_LEN + 1, b=np.infty, direction=lm.BACKWARD, g=g_bl, delta=SHAPE_LEN + 1,
                      label='segment 3')

# Defining the final cost function (a so called composite cost = CCost)
# mapping matrix between models and segments (rows = models, columns = segments)
F = [[0, 1, 0],
     [1, 1, 1]]
ccost = lm.CompositeCost((alssm_pulse, alssm_baseline), (segmentL, segmentC, segmentR), F)
print(ccost)

# filter signal
se_param = lm.SEParam(ccost, use_steady_state=False)
se_param.filter(y)  # run recursions

xs = se_param.minimize_x()  # unconstrained minimization

H_A = np.transpose(block_diag(np.eye(N2)[1:,:], np.eye(N1)))  # constrain matrix to find pulses of same shape as the reference pulse
H_0 = np.transpose(np.hstack([np.zeros((N1, N2)), np.eye(N1)]))  # constrain matrix to test for no pulse (baseline only)

print("H_A : ", H_A)
print("H_0 : ", H_0)

xs_A = se_param.minimize_x(H_A)
xs_0 = se_param.minimize_x(H_0)

J_A = se_param.eval_errors(xs_A, range(K))  # get SE (squared error) for hypothesis 1 (baseline + pulse)
J_0 = se_param.eval_errors(xs_0, range(K))  # get SE (squared error)  for hypothesis 0 (baseline only) --> J0 should be a vector not a matrice

LCR = -0.5 * np.log(J_A / J_0)  # log-cost ratio computation

vs = se_param.minimize_v(H_A)  # TODO: need to be joined with minimize_x
amp = vs[:, 0]

# find peaks
peaks, _ = find_peaks(LCR, height=LCR_THD, distance=MIN_DIST)


# --------------- plotting of results -----------------------
k = np.arange(K)

# Trajectories
trajs_baseline = lm.map_trajectories(ccost.trajectories(xs_A[peaks], F=[[0, 0, 0], [1, 1, 1]], thds=0.01), peaks, K,
                                     merge_ks=True, merge_seg=True)
trajs_pulse = lm.map_trajectories(ccost.trajectories(xs_A[peaks], F=[[0, 1, 0], [0, 1, 0]], thds=0.01), peaks, K,
                                  merge_ks=True, merge_seg=True)

fig, axs = plt.subplots(5, 1, sharex='all', figsize=(8, 6))

# Remove horizontal space between axes, maximize use of plotting pane
fig.tight_layout()
fig.subplots_adjust(hspace=0.0, left=0.08, bottom=0.05)

if peaks.size != 0:
    wins = lm.map_window(ccost.window(segment_selection=[True, True, True], thds=0.001), peaks, K, merge_ks=True,
                         fill_value=0)
    axs[0].plot(k, wins[0], color='r', lw=0.5, ls='-', label=segmentL.label)
    axs[0].plot(k, wins[1], color='g', lw=0.5, ls='-', label=segmentC.label)
    axs[0].plot(k, wins[2], color='b', lw=0.5, ls='-', label=segmentR.label)
axs[0].set(ylabel='windows')
axs[0].legend(loc='upper right')

axs[1].plot(k, y, color="gray", lw=1.0, label='y')
axs[1].legend(loc='upper right')

axs[2].plot(k, y, color="gray", lw=1.0, label='y')
axs[2].plot(range(K), trajs_baseline, color='g', lw=1.0, linestyle="-", label='traj. baseline')
axs[2].plot(range(K), trajs_pulse, color='r', lw=1.0, linestyle="-", label='traj. pulse')
axs[2].legend(loc='upper right')

axs[3].plot(k, J_A, lw=1.0, color='blue', label=r"$\tilde{J}(H v)$")
axs[3].plot(k, J_0, lw=1.0, color='darkred', label=r"$\tilde{J}(H_0 v)$")
axs[3].legend(loc='upper right')

axs[4].plot(k, LCR, lw=1.0, color='green', label=r"$LCR = -.5 ln(J(\hat{x}_A) / J(\hat{x}_0))$")
axs[4].plot(peaks, LCR[peaks], "x", color='r', markersize=8, markeredgewidth=2.0)
axs[4].axhline(LCR_THD, color="black", linestyle="--", lw=1.0)
axs[4].legend(loc='upper right')


plt.show()

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

Gallery generated by Sphinx-Gallery