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

Out:

(7, 4000)

import matplotlib.pyplot as plt
import numpy as np
import numpy.matlib
from scipy.signal import find_peaks
import lmlib as lm
from lmlib.utils.generator import gen_wgn, load_multi_channel

# --- Generating Test signal ---
K = 4000
NOISE = .2
fs = 600 # [Hz] sampling frequency

K_REF = 450 # reference index of shape to be detected along the signal

k = np.arange(K)

if True: # ecg test signal
    k_start = 68
    y1_mc = load_multi_channel('EECG_FILT_9CH_10S_FS600HZ.csv', K=K+k_start)
    y1 = y1_mc[k_start:,1] # select single channel
    yn = gen_wgn(K, .01, seed=1234)






# ----- Models and Parameter Estimation  --------------

# Polynomial ALSSM
alssm_poly = lm.AlssmPoly(poly_degree=3)
# Segment
segment = lm.Segment(a=-40, b=40, direction=lm.BACKWARD, g=5000)

# CostSeg
costs = lm.CostSegment(alssm_poly, segment)


# filter signal and take the approximation
y = y1 + yn*NOISE

se_param = lm.SEParam(costs)
se_param.filter(y)

xs = se_param.minimize_x() # unconstrained minimization
x_ref = xs[K_REF,:]  # get reference coefficients

H = np.reshape(x_ref, (-1,1))


vs = se_param.minimize_v(H) # unconstrained minimization

J_ref = se_param.eval_errors(vs@H.transpose())
J_0 = se_param.eval_errors(np.zeros_like(xs))

lcr = -.5*np.log(J_ref / J_0)

ind_peaks, _  = find_peaks(lcr, height=1.0, distance = fs*0.5, )

# ----------------  Plot  -----------------
ks = ind_peaks # indeces to display fit

trajs = lm.map_trajectories(costs.trajectories(xs[ks], thd=0.1), ks, K)
wins = lm.map_window(costs.window(), ks, K)


# axs[0].set( title='noise: $\sigma ='+str(NOISE)+'$')


fig = plt.figure(figsize=(7,5))

axs = fig.add_subplot(6, 1, 1)
print(wins.shape)
axs.plot(k, wins[0,:], lw=1, c='k', ls='--')
axs.set( ylabel='window(s)')


axs = fig.add_subplot(6, 1, (2,3), sharex=axs)
axs.plot(k, y, lw=1.0, c='k')
axs.set( ylabel='observations')

# for n in range(0,np.shape(trajs)[1]):
for (n, traj) in enumerate(trajs[:,:,0]):
    axs.plot(k, traj, lw=1.5, c='tab:green', ls='-', label='traj.' if n==0 else '')
#    axs.scatter(ks[n], traj[ks[n]], marker='x', c='tab:green')
for xc in ind_peaks:
    axs.axvline(x=xc, c='k', lw=0.5, ls='--')
axs.grid(True)
# axs[1].set_ylim([-3, 3])
axs.legend(loc='upper right')

axs = fig.add_subplot(6, 1, (4), sharex=axs)
axs.plot(k, J_ref, label='$J(Hv_k)$', c='tab:orange')
axs.plot(k, J_0, label='$J(0)$', c='tab:blue')
for xc in ind_peaks:
    axs.axvline(x=xc, c='k', lw=0.5, ls='--')
axs.legend(loc='upper right')
axs.set( ylabel='Costs')


axs = fig.add_subplot(6, 1, (5), sharex=axs)
axs.stem(ind_peaks, vs[ind_peaks], use_line_collection=True, label='$\lambda$')
for xc in ind_peaks:
    axs.axvline(x=xc, c='k', lw=0.5, ls='--')
axs.set_ylim([.8,1.2])
axs.axhline(y=1.0, c='k', lw=0.5, ls='--')
axs.set( ylabel='Amplitude $\lambda$')
axs.legend(loc='upper right')

axs = fig.add_subplot(6, 1, (6), sharex=axs)
axs.plot(k, lcr, label='LCR',  c='tab:purple')
axs.plot(ind_peaks, lcr[ind_peaks], "x", c='k')
for xc in ind_peaks:
    axs.axvline(x=xc, c='k', lw=0.5, ls='--')
axs.set( ylabel='LCR')

plt.subplots_adjust(hspace=0.4)

plt.show()

Peak detection with recursive moving average in ABP [ex9xx.x] Installers