Note
Click here to download the full example code
Out:
(4, 550)
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 = 550 NOISE = .05 OMEGA = 0.07 # frequency k = np.arange(K) if True: # ecg test signal k_start = 180 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 -------------- # Sinusoidal ALSSM alssm_poly = lm.AlssmSin(omega=OMEGA) # Segment segment = lm.Segment(a=0, b=int(np.floor(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 = [25, 140, 300, 500] # 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)) axs = fig.add_subplot(3, 1, 1) print(wins.shape) lss = ['--', '-', ':','-.','--'] for (n, _) in enumerate(ks): axs.plot(k, wins[n,:], lw=1, c='k', ls=lss[n]) axs.set( ylabel='window(s)') axs = fig.add_subplot(3, 1, (2,3)) axs.plot(k, y, lw=1.0, c='tab:gray') # for n in range(0,np.shape(trajs)[1]): for (n, traj) in enumerate(trajs[:,:,0]): axs.plot(k, traj, lw=1.5, c='k', ls='-', label='$a_0 cos(i \Omega) - a_1 sin(i \Omega)$' if n==0 else '') axs.scatter(ks[n], traj[ks[n]], marker='x', c='k') axs.plot(k, trajs0[n, :, 0], lw=1.5, c='k', ls='--', label='$a_0 cos(i \Omega)$' if n == 0 else '') axs.scatter(ks[n], trajs0[n, ks[n], 0], marker='x', c='k') axs.plot(k, trajs1[n, :, 0], lw=1.5, c='k', ls='-.', label='$-a_1 sin(i \Omega)$' if n == 0 else '') axs.scatter(ks[n], trajs1[n, ks[n], 0], marker='x', c='k') axs.grid(True) # axs[1].set_ylim([-3, 3]) axs.legend(loc='upper right') plt.subplots_adjust(hspace=0.4) plt.show()
Total running time of the script: ( 0 minutes 0.373 seconds)
Gallery generated by Sphinx-Gallery