Detection of changes in signal slope

Out:

[ 207  496  907 1108]

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, gen_slopes

# --- Generating Test signal ---


K = 2000
k = np.arange(K)



if True: # test signal with slope switches

    fs = 600  # [Hz] sampling frequency
    # create test signal
    y1 = gen_slopes(K, ks=[0, 200, 500, 900, 1100, 1700], deltas=[0, 0, 2.5, -0.5, -1, 0])
    yn = gen_wgn(K, sigma=0.05, seed=12345)

y = y1 + yn



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

# Polynomial ALSSM
alssm_line = lm.AlssmPoly(poly_degree=1)
# Segment

L=150
segment_left = lm.Segment(a=-L, b=-1, direction=lm.BACKWARD, g=5000)
segment_right = lm.Segment(a=0, b=L, direction=lm.BACKWARD, g=5000)
segment_left_right = lm.Segment(a=-L, b=L, direction=lm.BACKWARD, g=5000)


# CostSeg
cost_left = lm.CostSegment(alssm_line, segment_left)
cost_right = lm.CostSegment(alssm_line, segment_right)
cost_left_right = lm.CostSegment(alssm_line, segment_left_right)


# filter signal and take the approximation
se_param_left = lm.SEParam(cost_left)
se_param_left.filter(y)

se_param_right = lm.SEParam(cost_right)
se_param_right.filter(y)

se_param_left_right = lm.SEParam(cost_left_right)
se_param_left_right.filter(y)

# unconstrained minimizations
xs_left = se_param_left.minimize_x()
xs_right = se_param_right.minimize_x()
xs_left_right = se_param_left_right.minimize_x()

J_1_left = se_param_left.eval_errors(xs_left) # low if jump is detected
J_1_right = se_param_right.eval_errors(xs_right) # low if jump is detected
J_1 = J_1_left + J_1_right
J_0 = se_param_left_right.eval_errors(xs_left_right)  # low if no jumpm is detected

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

ind_peaks, _  = find_peaks(lcr, height=0.1, distance = 200, )

# ----------------  Plot  -----------------
ks = ind_peaks # indeces to display fit
print(ind_peaks)
trajs_left = lm.map_trajectories(cost_left.trajectories(xs_left[ks], thd=0.1), ks, K)
trajs_right = lm.map_trajectories(cost_right.trajectories(xs_right[ks], thd=0.1), ks, K)
trajs_left_right = lm.map_trajectories(cost_left_right.trajectories(xs_left_right[ks], thd=0.1), ks, K)



# fig = plt.figure(figsize=(7,3)) # for script
fig = plt.figure(figsize=(6,5)) # for presentation ETH


axs = fig.add_subplot(3, 1, (1))
axs.plot(k, y, lw=.5, c='tab:gray')
axs.plot(k, y1, ':', lw=1.0, c='k')
axs.set( ylabel='y')


for (n, _) in enumerate(trajs_left[:,:,0]):
    axs.plot(k, trajs_left[n,:,0], lw=1.5, c='b', ls='-', label='traj. $J_{a}^{-1}(\cdot)$' if n == 0 else '')
    axs.plot(k, trajs_right[n,:,0], lw=1.5, c='tab:green', ls='-', label='traj. $J_{0}^{b}(\cdot)$' if n == 0 else '')
    axs.plot(k, trajs_left_right[n, :, 0], lw=1.5, c='k', ls='--', label='traj. $J_{a}^{b}(\cdot)$' if n == 0 else '')
for xc in ind_peaks:
    axs.axvline(x=xc, c='b', lw=0.5, ls='--')
# axs.grid(True)
axs.legend(loc='upper right')

axs = fig.add_subplot(3, 1, (2), sharex=axs)
axs.plot(k, J_1_left, label='$J_{a}^{-1}(k,x_L)$', lw=1.0, ls='-', c='b')
axs.plot(k, J_1_right, label='$J_{0}^b(k,x_R)$', lw=1.0, ls='-', c='tab:green')
axs.plot(k, J_1, label='$H_1: J_{a}^{-1}(k,x_L)+J_{0}^b(k,x_R)$', lw=1.0, ls='-', c='k')
axs.plot(k, J_0, label='$H_0: J_{a}^{b}(k,x_L = x_R)$', lw=1.0, ls='--', c='k')

for xc in ind_peaks:
    axs.axvline(x=xc, c='b', lw=0.5, ls='--')
axs.legend(loc='upper right')
axs.set( ylabel='Costs')


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

axs.legend(loc='upper right')
plt.subplots_adjust(hspace=0.4)

plt.show()

Detection of changes in signal slope Peak detection with causal recursive moving average in ABP [ex9xx.x]