# Fitting of polynomial model to signal [ex901.0]ΒΆ

Out:

(3, 650)


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

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

sin_period = 800 # sinusoidal periodicity
k = np.arange(K)

if False: # sinusoidal test signal
y1 =  np.concatenate(
(
np.zeros(int(K/4)),
2*gen_sine(2*int(K/4), k_period=sin_period),
np.zeros(K-3*int(K/4))
))
yn = gen_wgn(K, 1, seed=1234)

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

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

# Sinusoidal ALSSM
alssm_poly = lm.AlssmPoly(poly_degree=2)
# Segment
segment = lm.Segment(a=0, b=70, 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)

H = [[1, 0, 0],
[0, 1, 0],
[0, 0, 1]]
xs = se_param.minimize_x(H) # unconstrained minimization

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

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

# ----------------  Plot  -----------------
ks = [50,  250,  450] # 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,3))

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

axs.plot(k, y, lw=1.0, c='k')

# 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='no constraint' if n==0 else '')
axs.scatter(ks[n], traj[ks[n]], marker='x', c='tab:green')
axs.plot(k, trajs0[n, :, 0], lw=1.5, c='tab:red', ls='-', label='1st order constr.' if n == 0 else '')
axs.scatter(ks[n], trajs0[n, ks[n], 0], marker='x', c='tab:red')
axs.plot(k, trajs1[n, :, 0], lw=1.5, c='tab:blue', ls='-', label='0 order constr.' if n == 0 else '')
axs.scatter(ks[n], trajs1[n, ks[n], 0], marker='x', c='tab:blue')
axs.grid(True)
# axs[1].set_ylim([-3, 3])
axs.legend(loc='upper right')