python - correlate - encontrar el cambio de tiempo entre dos formas de onda similares

cross correlation python

Question

Tengo que comparar dos formas de onda tiempo-voltaje. Debido a la peculiaridad de las fuentes de estas formas de onda, una de ellas puede ser una versión modificada en el tiempo de la otra.

¿Cómo puedo saber si hay un cambio de horario? Y si es así, ¿cuánto es?

Estoy haciendo esto en Python y deseo usar bibliotecas numpy / scipy.

Answer 1

Aquí hay otra opción:

from scipy import signal, fftpack

def get_max_correlation(original, match):
    z = signal.fftconvolve(original, match[::-1])
    lags = np.arange(z.size) - (match.size - 1)
    return ( lags[np.argmax(np.abs(z))] )

Answer 2

Depende del tipo de señal que tenga (¿periódica? ...), de si ambas señales tienen la misma amplitud y de qué precisión está buscando.

La función de correlación mencionada por highBandWidth podría funcionar para usted. Es lo suficientemente simple como para probarlo.

Otra opción más precisa es la que uso para el ajuste de líneas espectrales de alta precisión: usted modela su señal "maestra" con una spline y ajusta la señal de tiempo desplazado con la misma (aunque posiblemente sea necesario escalar la señal). Esto produce cambios de tiempo muy precisos. Una ventaja de este enfoque es que no tiene que estudiar la función de correlación. Por ejemplo, puedes crear la spline fácilmente con interpolate.UnivariateSpline() (desde SciPy). SciPy devuelve una función, que luego se ajusta fácilmente con optimize.leastsq ().

Answer 3

Esta función es probablemente más eficiente para señales de valor real. Utiliza rfft y cero rellena las entradas a una potencia de 2 lo suficientemente grande como para asegurar una correlación lineal (es decir, no circular):

def rfft_xcorr(x, y):
    M = len(x) + len(y) - 1
    N = 2 ** int(np.ceil(np.log2(M)))
    X = np.fft.rfft(x, N)
    Y = np.fft.rfft(y, N)
    cxy = np.fft.irfft(X * np.conj(Y))
    cxy = np.hstack((cxy[:len(x)], cxy[N-len(y)+1:]))
    return cxy

El valor de retorno es la longitud M = len(x) + len(y) - 1 (hackeado junto con hstack para eliminar los ceros adicionales del redondeo hasta una potencia de 2). Los retardos no negativos son cxy[0], cxy[1], ..., cxy[len(x)-1] , mientras que los retardos negativos son cxy[-1], cxy[-2], ..., cxy[-len(y)+1] .

Para hacer coincidir una señal de referencia, rfft_xcorr(x, ref) y busco el pico. Por ejemplo:

def match(x, ref):
    cxy = rfft_xcorr(x, ref)
    index = np.argmax(cxy)
    if index < len(x):
        return index
    else: # negative lag
        return index - len(cxy)   

In [1]: ref = np.array([1,2,3,4,5])
In [2]: x = np.hstack(([2,-3,9], 1.5 * ref, [0,3,8]))
In [3]: match(x, ref)
Out[3]: 3
In [4]: x = np.hstack((1.5 * ref, [0,3,8], [2,-3,-9]))
In [5]: match(x, ref)
Out[5]: 0
In [6]: x = np.hstack((1.5 * ref[1:], [0,3,8], [2,-3,-9,1]))
In [7]: match(x, ref)
Out[7]: -1

No es una forma robusta de hacer coincidir las señales, pero es rápida y fácil.

Answer 4

Si uno está desplazado en el tiempo por el otro, verá un pico en la correlación. Como el cálculo de la correlación es costoso, es mejor usar FFT. Entonces, algo como esto debería funcionar:

af = scipy.fft(a)
bf = scipy.fft(b)
c = scipy.ifft(af * scipy.conj(bf))

time_shift = argmax(abs(c))

Answer 5

scipy proporciona una función de correlación que funcionará bien para una entrada pequeña y también si desea una correlación no circular, lo que significa que la señal no se ajustará. tenga en cuenta que en mode=''full'' , el tamaño de la matriz devuelta por signal.correlation es la suma de los tamaños de señal de entrada - 1, por lo que el valor de argmax está desactivado por (tamaño de señal -1 = 20) de lo que parece esperar.

from scipy import signal, fftpack
import numpy
a = numpy.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0])
b = numpy.array([0, 0, 0, 0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0])
numpy.argmax(signal.correlate(a,b)) -> 16
numpy.argmax(signal.correlate(b,a)) -> 24

Los dos valores diferentes corresponden a si el desplazamiento está en a o b .

Si desea una correlación circular y para un gran tamaño de señal, puede usar el teorema de convolución / transformada de Fourier con la advertencia de que la correlación es muy similar pero no idéntica a la de convolución.

A = fftpack.fft(a)
B = fftpack.fft(b)
Ar = -A.conjugate()
Br = -B.conjugate()
numpy.argmax(numpy.abs(fftpack.ifft(Ar*B))) -> 4
numpy.argmax(numpy.abs(fftpack.ifft(A*Br))) -> 17

de nuevo, los dos valores corresponden a si su interpretación de un cambio en a o un cambio en b .

La conjugación negativa se debe a la conversión por convolución de una de las funciones, pero en la correlación no hay ningún cambio. Puede deshacer el cambio de posición invirtiendo una de las señales y luego tomando la FFT, o tomando la FFT de la señal y luego tomando el conjugado negativo. es decir, lo siguiente es cierto: Ar = -A.conjugate() = fft(a[::-1])