Come trovare zero incroci con isteresi?

Question

In numpy, vorrei rilevare i punti dai quali il segnale attraversa (essendo stato precedentemente) al di sotto di una certa soglia, per essere sopra una certa altra soglia. Questo è per cose come debouncing, o accurati zero crossing in presenza di rumore, ecc

Ti piace questa:

import numpy 

# set up little test problem 
N = 1000 
values = numpy.sin(numpy.linspace(0, 20, N)) 
values += 0.4 * numpy.random.random(N) - 0.2 
v_high = 0.3 
v_low = -0.3 

# find transitions from below v_low to above v_high  
transitions = numpy.zeros_like(values, dtype=numpy.bool) 

state = "high" 

for i in range(N): 
    if values[i] > v_high: 
     # previous state was low, this is a low-to-high transition 
     if state == "low": 
      transitions[i] = True 
     state = "high" 
    if values[i] < v_low: 
     state = "low" 

vorrei un modo per fare questo senza loop sopra la matrice esplicitamente: ma io non posso pensare in alcun modo, dal momento che ogni valore di stato dipende dallo stato precedente. È possibile fare a meno del ciclo?

Answer 1

Questo può essere fatto in questo modo:

def hyst(x, th_lo, th_hi, initial = False): 
    hi = x >= th_hi 
    lo_or_hi = (x <= th_lo) | hi 
    ind = np.nonzero(lo_or_hi)[0] 
    if not ind.size: # prevent index error if ind is empty 
     return np.zeros_like(x, dtype=bool) | initial 
    cnt = np.cumsum(lo_or_hi) # from 0 to len(x) 
    return np.where(cnt, hi[ind[cnt-1]], initial) 

Spiegazione: ind sono gli indici di tutti i campioni in cui il segnale è al di sotto del più basso o sopra la soglia superiore, e per cui la posizione del 'interruttore 'è quindi ben definito. Con cumsum, si crea una sorta di contatore che punta all'indice dell'ultimo campione ben definito. Se l'inizio del vettore di input è compreso tra le due soglie, cnt sarà 0, quindi è necessario impostare l'uscita corrispondente sul valore iniziale utilizzando la funzione where.

Credito: questo è un trucco che ho trovato in un old post su un forum Matlab, che ho tradotto in Numpy. Questo codice è un po 'difficile da capire e ha anche bisogno di allocare vari array intermedi. Sarebbe meglio se Numpy includesse una funzione dedicata, simile al proprio ciclo for, ma implementata in C per la velocità.

prova rapida:

x = np.linspace(0,20, 1000) 
y = np.sin(x) 
h1 = hyst(y, -0.5, 0.5) 
h2 = hyst(y, -0.5, 0.5, True) 
plt.plot(x, y, x, -0.5 + h1, x, -0.5 + h2) 
plt.legend(('input', 'output, start=0', 'output, start=1')) 
plt.title('Thresholding with hysteresis') 
plt.show() 

Risultato:

Answer 2

Le modifiche che ho dovuto fare per il mio lavoro, tutti basati sulla risposta di cui sopra da Bas Swinckels, per consentire il rilevamento di soglia di passaggio per quando si utilizza di serie nonché soglie invertite.

io non sono felice con la denominazione duro, forse ora dovrebbe leggere th_hi2lo e th_lo2hi invece di th_lo e th_hi? Usando i valori originali, il comportamento è lo stesso duro.

def hyst(x, th_lo, th_hi, initial = False): 
    """ 
    x : Numpy Array 
     Series to apply hysteresis to. 
    th_lo : float or int 
     Below this threshold the value of hyst will be False (0). 
    th_hi : float or int 
     Above this threshold the value of hyst will be True (1). 
    """   

    if th_lo > th_hi: # If thresholds are reversed, x must be reversed as well 
     x = x[::-1] 
     th_lo, th_hi = th_hi, th_lo 
     rev = True 
    else: 
     rev = False 

    hi = x >= th_hi 
    lo_or_hi = (x <= th_lo) | hi 

    ind = np.nonzero(lo_or_hi)[0] # Index für alle darunter oder darüber 
    if not ind.size: # prevent index error if ind is empty 
     x_hyst = np.zeros_like(x, dtype=bool) | initial 
    else: 
     cnt = np.cumsum(lo_or_hi) # from 0 to len(x) 
     x_hyst = np.where(cnt, hi[ind[cnt-1]], initial) 

    if rev: 
     x_hyst = x_hyst[::-1] 

    return x_hyst 

e come sopra un test del codice per vedere che cosa fa:

x = np.linspace(0,20, 1000) 
y = np.sin(x) 
h1 = hyst(y, -0.2, 0.2) 
h2 = hyst(y, +0.5, -0.5) 
plt.plot(x, y, x, -0.2 + h1*0.4, x, -0.5 + h2) 
plt.legend(('input', 'output, classic, hyst(y, -0.2, +0.2)', 
      'output, reversed, hyst(y, +0.5, -0.5)')) 
plt.title('Thresholding with hysteresis') 
plt.show()