Estoy usando cython para un cálculo de correlación en mi programa python. Tengo dos conjuntos de datos de audio y necesito saber la diferencia horaria entre ellos. El segundo conjunto se corta según los tiempos de inicio y luego se desliza por el primer conjunto. Hay dos bucles for: uno desliza el conjunto y el bucle interno calcula la correlación en ese punto. Este método funciona muy bien y es lo suficientemente preciso.

El problema es que con Python puro esto lleva más de un minuto. Con mi código de cython, lleva unos 17 segundos. Esto todavía es demasiado. ¿Tiene alguna pista sobre cómo acelerar este código:

import numpy as np
cimport numpy as np

cimport cython

FTYPE = np.float
ctypedef np.float_t FTYPE_t

@cython.boundscheck(False)
def delay(np.ndarray[FTYPE_t, ndim=1] f, np.ndarray[FTYPE_t, ndim=1] g):
    cdef int size1 = f.shape[0]
    cdef int size2 = g.shape[0]
    cdef int max_correlation = 0
    cdef int delay = 0
    cdef int current_correlation, i, j

    # Move second data set frame by frame
    for i in range(0, size1 - size2):
        current_correlation = 0

        # Calculate correlation at that point
        for j in range(size2):
            current_correlation += f[<unsigned int>(i+j)] * g[j]

        # Check if current correlation is highest so far
        if current_correlation > max_correlation:
            max_correlation = current_correlation
            delay = i

    return delay
16
jushie 29 jul. 2009 a las 16:39

3 respuestas

La mejor respuesta

Editar:
Ahora hay scipy.signal.fftconvolve que sería el enfoque preferido para hacer el enfoque de convolución basado en FFT que describo a continuación. Dejaré la respuesta original para explicar el problema de la velocidad, pero en la práctica use scipy.signal.fftconvolve.

Respuesta original:
Utilizando FFTs y el teorema de convolución le dará ganancias de velocidad dramáticas al convertir el problema de O (n ^ 2) a O (n log n). Esto es particularmente útil para conjuntos de datos largos, como el suyo, y puede proporcionar ganancias de velocidad de 1000 o mucho más, dependiendo de la longitud. También es fácil de hacer: solo FFT ambas señales, multiplicar e invertir FFT el producto. numpy.correlate no usa el método FFT en la rutina de correlación cruzada y se usa mejor con núcleos muy pequeños.

Aquí hay un ejemplo

from timeit import Timer
from numpy import *

times = arange(0, 100, .001)

xdata = 1.*sin(2*pi*1.*times) + .5*sin(2*pi*1.1*times + 1.)
ydata = .5*sin(2*pi*1.1*times)

def xcorr(x, y):
    return correlate(x, y, mode='same')

def fftxcorr(x, y):
    fx, fy = fft.fft(x), fft.fft(y[::-1])
    fxfy = fx*fy
    xy = fft.ifft(fxfy)
    return xy

if __name__ == "__main__":
    N = 10
    t = Timer("xcorr(xdata, ydata)", "from __main__ import xcorr, xdata, ydata")
    print 'xcorr', t.timeit(number=N)/N
    t = Timer("fftxcorr(xdata, ydata)", "from __main__ import fftxcorr, xdata, ydata")
    print 'fftxcorr', t.timeit(number=N)/N

Lo que proporciona los tiempos de ejecución por ciclo (en segundos, para una forma de onda larga de 10.000)

xcorr 34.3761689901
fftxcorr 0.0768054962158

Está claro que el método fftxcorr es mucho más rápido.

Si traza los resultados, verá que son muy similares cerca del cambio de tiempo cero. Sin embargo, tenga en cuenta que a medida que se aleje, el xcorr disminuirá y el fftxcorr no lo hará. Esto se debe a que es un poco ambiguo qué hacer con las partes de la forma de onda que no se superponen cuando las formas de onda se desplazan. xcorr lo trata como cero y el FFT trata las formas de onda como periódicas, pero si se trata de un problema, se puede solucionar con relleno de cero.

37
tom10 12 oct. 2014 a las 16:29
  • puede extraer el rango (size2) del bucle externo
  • puede usar sum () en lugar de un bucle para calcular la correlación actual
  • puede almacenar correlaciones y demoras en una lista y luego usar max () para obtener la más grande
2
dugres 29 jul. 2009 a las 13:04

El truco con este tipo de cosas es encontrar una manera de dividir y conquistar.

Actualmente, se está deslizando a cada posición y verifica cada punto en cada posición, efectivamente una operación O (n ^ 2).

Debe reducir la comprobación de cada punto y la comparación de cada posición a algo que trabaje menos para determinar una no coincidencia.

Por ejemplo, podría tener un breve "¿está esto incluso cerca?" filtro que verifica las primeras posiciones. Si la correlación está por encima de algún umbral, entonces continúe, de lo contrario, renuncie y siga adelante.

Podría tener un "cheque cada 8a posición" que multiplique por 8. Si esto es demasiado bajo, omítalo y continúe. Si esto es lo suficientemente alto, verifique todos los valores para ver si ha encontrado los máximos.

El problema es el tiempo requerido para hacer todas estas multiplicaciones: (f[<unsigned int>(i+j)] * g[j]) En efecto, está llenando una matriz grande con todos estos productos y seleccionando la fila con la suma máxima. No desea calcular "todos" los productos. Solo lo suficiente de los productos para asegurarse de haber encontrado la suma máxima.

El problema con encontrar máximos es que tienes que sumar todo para ver si es más grande. Si puede convertir esto en un problema de minimización, es más fácil abandonar los productos informáticos y las sumas una vez que un resultado intermedio supera un umbral.

(Creo que esto podría funcionar. No lo he probado).

Si utilizó max(g)-g[j] para trabajar con números negativos, estaría buscando el más pequeño, no el más grande. Podrías calcular la correlación para la primera posición. Cualquier cosa que se sumara a un valor mayor podría detenerse de inmediato: no más multiplicaciones o sumas para ese desplazamiento, cambie a otro.

2
S.Lott 29 jul. 2009 a las 12:58