# Constrained fit of sinusoidal models to ECG signal [ex901.0]ΒΆ

Out:

(3, 600)

import matplotlib.pyplot as plt
import numpy as np
import lmlib as lm
from lmlib.utils.generator import gen_sine, gen_rectangle, gen_wgn, load_multi_channel, k_period_to_omega

# --- Generating Test signal ---
K = 600
NOISE = .05

OMEGA = 0.05 # sine wave frequency
k = np.arange(K)

if True: # sinusoidal test signal
y1 =  np.concatenate(
(
np.zeros(int(K/8)),
-1+2*gen_rectangle(3*int(K/4), k_period=125, k_on = 62),
np.zeros(K-7*int(K/8))
))
yn = gen_wgn(K, 1, seed=1234)

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

# Sinusoidal ALSSM
alssm_poly = lm.AlssmSin(omega=OMEGA)
# Segment
segment = lm.Segment(a=0, b=int(np.floor(2*np.pi/OMEGA)), 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

H = [[0],
[1]]
xs0 = se_param.minimize_x(H) # constrained minimization

H = [[1],
[0]]
xs1 = se_param.minimize_x(H) # constrained minimization

# ----------------  Plot  -----------------
ks = [50, 200,  350, ] # indeces to display fit

trajs = lm.map_trajectories(costs.trajectories(xs[ks], thd=0.1), ks, K)
trajs0 = lm.map_trajectories(costs.trajectories(xs0[ks], thd=0.1), ks, K)
trajs1 = lm.map_trajectories(costs.trajectories(xs1[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=(8,6))

print(wins.shape)
lss = ['-', '--', ':','-.','--']
for (n, _) in enumerate(ks):
axs.plot(k, wins[n,:], lw=1, c='tab:gray', ls=lss[n])

axs.set( ylabel='window(s)')

for (n, traj) in enumerate(trajs[:,:,0]):
axs.axhline(y=0, lw=.75, ls='--', c='k')
axs.plot(k, y, lw=.75, c='tab:gray')
axs.plot(k, traj, lw=1.5, c='k', ls='-', label=r'$v \cdot sin(i \Omega + \Phi)$' if n==0 else '')
axs.scatter(ks[n], traj[ks[n]], marker='x', c='k')
axs.legend(loc='upper right')

for (n, traj) in enumerate(trajs0[:, :, 0]):
axs.axhline(y=0, lw=.75, ls='--', c='k')
axs.plot(k, y, lw=.75, c='tab:gray')
axs.plot(k, traj, lw=1.5, c='k', ls='-', label=r'$v \cdot sin(i \Omega + 0)$' if n == 0 else '')
axs.scatter(ks[n], traj[ks[n]], marker='x', c='k')
axs.legend(loc='upper right')

for (n, traj) in enumerate(trajs1[:, :, 0]):
axs.axhline(y=0, lw=.75, ls='--', c='k')
axs.plot(k, y, lw=.75, c='tab:gray')
axs.plot(k, traj, lw=1.5, c='k', ls='-', label=r'$v \cdot sin(i \Omega + pi/2)$' if n == 0 else '')
axs.scatter(ks[n], traj[ks[n]], marker='x', c='k')
axs.legend(loc='upper right')