Detection of jumps in a signal [ex911.0]ΒΆ

../../_images/sphx_glr_fig-examp-jumpdetector-noise_001.png
import numpy as np

import matplotlib.pyplot as plt
import lmlib as lm

from scipy.signal import find_peaks
from lmlib.utils.generator import *

from L1_costs import estAverageRect # import cost computation functions
from L1_costs import calcCostRatio
from L1_costs import getSEParams
from L1_costs import findPeaksLCR

from L1_plots import plotFigsCR3 # import plotting functions




seed = 1234578 # set pseudo-radom generator to fixed start value --> leads to reproducible results

NOF_X = 100
J_MAX = 100
a = -100 # rectangular window
b = 100
K = 2000 # number of samples
k = np.arange(0, K)
gamma = .992 # decay factor
min_distance = 150 # minimal distance between two detected peaks
min_lcr = 0.025 # minimal lcr to detect an event

NOISES = [.05,0.25,1.5] # Displaying multiple noise levels


# generate windows
i_w = np.arange(-K,K) # indices of window weights
w_Left = np.zeros(2*K); # generate window
np.putmask(w_Left, i_w<=0, gamma**np.abs(i_w))

w_Right = np.zeros(2*K); # generate window
np.putmask(w_Right, i_w>0, gamma**np.abs(i_w))



# create test signal
y1 = gen_steps(K, [500, 700, 1100, 1400], [1, .5, -1, .2])
yn = gen_wgn(K, sigma=1, seed=seed)

y = y1+yn


# ------------
# Jump detection on examples of different noise levels
y_s = []
CR_s = []
i_peaks_s = []
dx_s = []
labels = []
# ii = 10
for noise in NOISES:
    y = y1 + yn*noise
    (W_Left, xi_Left, kappa_Left) = getSEParams(y, a,b, w_Left, i_w)
    (W_Right, xi_Right, kappa_Right) = getSEParams(y, a,b, w_Right, i_w)

    (CR,x_hatLeft, x_hatRight, x_hat, J1, J2, J12) = calcCostRatio(W_Left, xi_Left, kappa_Left, W_Right, xi_Right, kappa_Right)

    dx = x_hatRight-x_hatLeft
    (i_peaks, dx_peaks) = findPeaksLCR(-.5*np.log(CR), dx, min_distance = min_distance, min_val = min_lcr)
    y_s.append(y)
    CR_s.append(CR)
    i_peaks_s.append(i_peaks)
    dx_s.append(dx)
    labels.append("noise="+str(noise))

# multiple plots for multiple nois levels
plotFigsCR3(y_s, k, CR_s, i_peaks_s, dx_s, labels)

plt.show()

Total running time of the script: ( 0 minutes 6.209 seconds)

Gallery generated by Sphinx-Gallery