Correlazione incrociata della funzione non periodica con NumPy

Question

Ho due serie di dati che sto cercando di correlare. Sembrano simili alla funzione arctan, quindi l'ho usata come modello per capire come eseguire l'elaborazione del segnale.

x = linspace(-15, 15, 2**13) 
f1 = arctan(x) 
f2 = arctan(x + 2) 

La domanda che ho bisogno di rispondere è, quanto ho bisogno di spostare il segnale verde per farlo (per lo più) sovrapposizione con quello blu? Ho pensato che sarebbe stato semplice trovare il massimo nella funzione di cross-correlazione di f1 e f2, e ho seguito ampiamente il consiglio qui: How to correlate two time series with gaps and different time bases?. Questo è quello che ho cercato

c = correlate(f1, f2, 'full') 
s = arange(1-2**13, 2**13) 
dx = 30/2**13 
shift = s[c.argmax()]*dx 

mi aspetterei shift di eguagliare più o meno esattamente 2, ma in realtà è solo 0.234. Questo non ha alcun senso per me; Ho trovato la coordinata x del massimo nella correlazione incrociata, che dovrebbe essere trovata dove c'è la massima sovrapposizione dei due segnali.

Qualche idea su come calcolare questa quantità per questo tipo di funzione?

EDIT: Vorrei aggiungere, per i miei dati reali, tutti i valori saranno tra zero e uno

EDIT EDIT: Le seguenti funzioni sono in realtà più come il mio vero dati:

x = linspace(-15, 15, 400) 
f1 = (arctan(-x) + pi/2)/pi 
f2 = (arctan(-x + 2) + pi/2)/pi 

Quindi, utilizzando la formula qui data: http://paulbourke.net/miscellaneous/correlate/ posso scrivere una funzione di correlazione incrociata che pads i dati da aggiungere quelli a sinistra e zeri a destra:

def xcorr(x, y); 
    mx = x.mean() 
    my = y.mean() 
    sx = x.std() 
    sy = y.std() 
    r = zeros(2*len(x)) 

    for d in range(-len(x), len(x)): 
     csum = 0 
     for i in range(0, len(x): 
      yindx = i - d 
      if i - d < 0: 
       yval = 1 
      elif i - d >= len(x): 
       yval = 0 
      else: 
       yval = y[yindx] 
      csum += (x[i] - mx) * (yval - my) 
     r[d + len(x)] = csum/(sx * sy) 
    return r 

Con questa funzione, ora posso fare

c = xcorr(f1, f2) 
s = arange(-400, 400) 
dx = 30/400 
shift = s[c.argmax()] * dx 

che esce a 2.025, che è il più vicino si può arrivare a 2 con questa precisione. Quindi sembra che Jamie avesse ragione, il problema è in che modo il numding di correlate fa il riempimento dei segnali.

Quindi, ovviamente, la mia funzione xcorr è molto lenta così com'è. La domanda ora è, c'è un modo per fare in modo che NumPy faccia una cosa simile, o dovrei semplicemente scrivere il mio algoritmo in C usando ctypes?

Answer 1

Ciò che il massimo della correlazione incrociata rileva è lo spostamento a cui la somma dei prodotti dei due segnali è massima. Si potrebbe pensare che la correlazione incrociata di due segnali, uno uno nel tempo dell'altra, sarebbe massima al turno. E mentre questo è vero per i segnali infiniti, numpy zero-pad i tuoi due segnali, quindi con due segnali ugualmente lunghi, hai solo una sommatoria su 2**13 valori diversi da zero per uno spostamento di zero, e i valori più alti da una migliore corrispondenza di le due funzioni, sono sfalsate ma ci sono meno valori diversi da zero.

Se i segnali erano 0 a +/- infinito, questo non sarebbe male. Così com'è, non sono stato in grado di trovare alcuna soluzione ragionevole usando la cross-correlazione.

Answer 2

Come le persone hanno sottolineato, la correlazione incrociata è confusa dal riempimento oltre i dati.

Anche se sembra che tu stia buttando via dati buoni, spesso è meglio tagliare semplicemente il set di dati in modo che la correlazione possa essere fatta senza supposizioni (almeno se confrontato con l'alternativa di correlare i dati effettivi con i dati composti per l'imbottitura).

x = linspace(-15, 15, 4000) 
f1 = (arctan(-x) + pi/2)/pi 
f2 = (arctan(-x + 2) + pi/2)/pi 

L4 = int(len(f2)/8) 
sf2 = f2[L4:-L4] 

c = correlate(f1-mean(f1), sf2-mean(f1), 'same') 
print "peak correlation occurs at:", x[argmax(c)] # -2.02925731433 

plt.plot(x, c) 
plt.show() 

Inoltre, però, non sono sicuro che xcorr è l'approccio migliore qui. Che ne dici di sommare la distanza tra i valori dell'asse y per diversi turni e prendere il minimo, che si allontana da tutti i problemi di dove lo zero è, ecc.