¿Cómo limitar el ancho de la ventana de correlación cruzada en Numpy?

Estoy aprendiendo entumecimiento / scipy, proveniente de un fondo de MATLAB. La función xcorr en Matlab tiene un argumento opcional “maxlag” que limita el intervalo de demora de –maxlag a maxlag. Esto es muy útil si está observando la correlación cruzada entre dos series de tiempo muy largas pero solo está interesado en la correlación dentro de un cierto rango de tiempo. Los aumentos de rendimiento son enormes, considerando que la correlación cruzada es increíblemente costosa de calcular.

En numpy / scipy parece que hay varias opciones para calcular la correlación cruzada. numpy.correlate , numpy.convolve , scipy.signal.fftconvolve . Si alguien desea explicar la diferencia entre estos, me encantaría escucharlo, pero principalmente lo que me preocupa es que ninguno de ellos tiene una función de máximo rendimiento. Esto significa que incluso si solo quiero ver correlaciones entre dos series de tiempo con retrasos entre -100 y +100 ms, por ejemplo, todavía calculará la correlación para cada retraso entre -20000 y +20000 ms (que es la longitud de la serie de tiempo). Esto da un golpe de rendimiento de 200x! ¿Tengo que recodificar la función de correlación cruzada a mano para incluir esta característica?

Aquí hay un par de funciones para calcular la correlación automática y cruzada con retrasos limitados. El orden de multiplicación (y conjugación, en el caso complejo) se eligió para coincidir con el comportamiento correspondiente de numpy.correlate .

 import numpy as np from numpy.lib.stride_tricks import as_strided def _check_arg(x, xname): x = np.asarray(x) if x.ndim != 1: raise ValueError('%s must be one-dimensional.' % xname) return x def autocorrelation(x, maxlag): """ Autocorrelation with a maximum number of lags. `x` must be a one-dimensional numpy array. This computes the same result as numpy.correlate(x, x, mode='full')[len(x)-1:len(x)+maxlag] The return value has length maxlag + 1. """ x = _check_arg(x, 'x') p = np.pad(x.conj(), maxlag, mode='constant') T = as_strided(p[maxlag:], shape=(maxlag+1, len(x) + maxlag), strides=(-p.strides[0], p.strides[0])) return T.dot(p[maxlag:].conj()) def crosscorrelation(x, y, maxlag): """ Cross correlation with a maximum number of lags. `x` and `y` must be one-dimensional numpy arrays with the same length. This computes the same result as numpy.correlate(x, y, mode='full')[len(a)-maxlag-1:len(a)+maxlag] The return vaue has length 2*maxlag + 1. """ x = _check_arg(x, 'x') y = _check_arg(y, 'y') py = np.pad(y.conj(), 2*maxlag, mode='constant') T = as_strided(py[2*maxlag:], shape=(2*maxlag+1, len(y) + 2*maxlag), strides=(-py.strides[0], py.strides[0])) px = np.pad(x, maxlag, mode='constant') return T.dot(px) 

Por ejemplo,

 In [367]: x = np.array([2, 1.5, 0, 0, -1, 3, 2, -0.5]) In [368]: autocorrelation(x, 3) Out[368]: array([ 20.5, 5. , -3.5, -1. ]) In [369]: np.correlate(x, x, mode='full')[7:11] Out[369]: array([ 20.5, 5. , -3.5, -1. ]) In [370]: y = np.arange(8) In [371]: crosscorrelation(x, y, 3) Out[371]: array([ 5. , 23.5, 32. , 21. , 16. , 12.5, 9. ]) In [372]: np.correlate(x, y, mode='full')[4:11] Out[372]: array([ 5. , 23.5, 32. , 21. , 16. , 12.5, 9. ]) 

(Será bueno tener una característica de este tipo en el propio numpy).

matplotlib.pyplot proporciona una syntax similar a matlab para calcular y trazar la correlación cruzada, la correlación automática, etc.

Puede usar xcorr que permite definir el parámetro maxlags .

  import matplotlib.pyplot as plt import numpy as np data = np.arange(0,2*np.pi,0.01) y1 = np.sin(data) y2 = np.cos(data) coeff = plt.xcorr(y1,y2,maxlags=10) print(*coeff) [-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10] [ -9.81991753e-02 -8.85505028e-02 -7.88613080e-02 -6.91325329e-02 -5.93651264e-02 -4.95600447e-02 -3.97182508e-02 -2.98407146e-02 -1.99284126e-02 -9.98232812e-03 -3.45104289e-06 9.98555430e-03 1.99417667e-02 2.98641953e-02 3.97518558e-02 4.96037706e-02 5.94189688e-02 6.91964864e-02 7.89353663e-02 8.86346584e-02 9.82934198e-02]  Line2D(_line0) 

Hasta que numpy implemente el argumento maxlag, puede usar la función ucorrelate del paquete pycorrelate . ucorrelate opera en matrices numpy y tiene una palabra clave maxlag . Implementa la correlación usando un for-loop y optimiza la velocidad de ejecución con numba.

Ejemplo : autocorrelación con 3 intervalos de tiempo:

 import numpy as np import pycorrelate as pyc x = np.array([2, 1.5, 0, 0, -1, 3, 2, -0.5]) c = pyc.ucorrelate(x, x, maxlag=3) c 

Resultado:

 Out[1]: array([20, 5, -3]) 

La documentación de pycorrelate contiene un cuaderno que muestra una coincidencia perfecta entre pycorrelate.ucorrelate y numpy.correlate :

introduzca la descripción de la imagen aquí

Creo que he encontrado una solución, ya que me enfrentaba al mismo problema:

Si tienes dos vectores y de cualquier longitud N, y quieres una correlación cruzada con una ventana de len m fija, puedes hacer:

 x =  y =  # Trim your variables x_short = x[window:] y_short = y[window:] # do two xcorrelations, lagging x and y respectively left_xcorr = np.correlate(x, y_short) #defaults to 'valid' right_xcorr = np.correlate(x_short, y) #defaults to 'valid' # combine the xcorrelations # note the first value of right_xcorr is the same as the last of left_xcorr xcorr = np.concatenate(left_xcorr, right_xcorr[1:]) 

Recuerde que es posible que necesite normalizar las variables si desea una correlación acotada

Aquí hay otra respuesta, obtenida de aquí , parece más rápida en el margen que np.correlate y tiene el beneficio de devolver una correlación normalizada:

 def rolling_window(self, a, window): shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) strides = a.strides + (a.strides[-1],) return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) def xcorr(self, x,y): N=len(x) M=len(y) meany=np.mean(y) stdy=np.std(np.asarray(y)) tmp=self.rolling_window(np.asarray(x),M) c=np.sum((y-meany)*(tmp-np.reshape(np.mean(tmp,-1),(N-M+1,1))),-1)/(M*np.std(tmp,-1)*stdy) return c 

Como respondí aquí, https://stackoverflow.com/a/47897581/5122657 matplotlib.xcorr tiene el parámetro maxlags. En realidad, es una envoltura de numpy.correlate , por lo que no hay ahorro de rendimiento. Sin embargo, da exactamente el mismo resultado dado por la función de correlación cruzada de Matlab. A continuación edité el código de matplotlib para que devuelva solo la correlación. La razón es que si usamos matplotlib.corr como es, también devolverá la ttwig. El problema es que si ponemos un tipo de datos complejo como los argumentos en él, obtendremos una advertencia de “conversión compleja a tipo de datos real” cuando matplotlib intente dibujar la ttwig.

  import numpy as np import matplotlib.pyplot as plt def xcorr(x, y, maxlags=10): Nx = len(x) if Nx != len(y): raise ValueError('x and y must be equal length') c = np.correlate(x, y, mode=2) if maxlags is None: maxlags = Nx - 1 if maxlags >= Nx or maxlags < 1: raise ValueError('maxlags must be None or strictly positive < %d' % Nx) c = c[Nx - 1 - maxlags:Nx + maxlags] return c