Cython e la velocità numpy

Question

sto usando Cython per un calcolo di correlazione nel mio programma python. Ho due insiemi di dati audio e ho bisogno di sapere la differenza di tempo tra di loro. Il secondo set è tagliato in base a tempi di inizio e poi scivolò sul primo set. Ci sono due cicli for: uno scorre il set e il ciclo interno calcola correlazione in quel punto. Questo metodo funziona molto bene ed è abbastanza preciso.

Il problema è che con puro python questo richiede più di un minuto. Con il mio codice Cython, ci vogliono circa 17 secondi. Questo è ancora troppo. Avete suggerimenti su come accelerare questo codice:

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

Answer 1

Modifica
C'è ora scipy.signal.fftconvolve che sarebbe il preferito approccio al fare l'approccio basato convoluzione FFT che descrivo qui di seguito. Lascio la risposta originale per spiegare la questione di velocità, ma in uso la pratica scipy.signal.fftconvolve.

risposta originale:
Uso di FFT e circonvoluzione teorema vi darà guadagni di velocità drammatici convertendo il problema da O (n ^ 2) a O (n log n). Ciò è particolarmente utile per insiemi di dati lunghi, come la vostra, e può dare guadagni di velocità di 1000 o più, a seconda della lunghezza. E 'anche facile da fare: basta FFT entrambi i segnali, si moltiplicano, e FFT inversa del prodotto. numpy.correlate non utilizza il metodo FFT nella routine correlazione incrociata ed è meglio utilizzato con molto piccole kernel.

Ecco un esempio

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

che dà i tempi di esecuzione per ciclo (in secondi, per un lungo forma d'onda di 10.000)

xcorr 34.3761689901
fftxcorr 0.0768054962158

E 'chiaro il metodo fftxcorr è molto più veloce.

Se si traccia i risultati, vedrai che sono molto simili nei pressi di spostamento di tempo zero. Si noti, però, come si ottiene più lontano il xcorr diminuisce e la fftxcorr non lo farà. Questo è perché è un po 'ambigua cosa fare con le parti della forma d'onda che non si sovrappongono quando le forme d'onda sono spostati. xcorr tratta come zero e la FFT tratta le forme d'onda periodica, ma se è un problema può essere risolto da zero padding.

Answer 2

Il trucco con questo genere di cose è quello di trovare un modo per dividere e conquistare.

Al momento, si sta scorrevole per ogni posizione e controllare ogni punto in ogni posizione -. Efficacemente un O (n ^ 2) il funzionamento

È necessario ridurre il controllo di tutti punto e il confronto delle tutti grado di qualcosa che fa meno lavoro per determinare una non-partita.

Per esempio, si potrebbe avere una più breve "è presente anche vicino?" filtro che controlla le prime posizioni. Se la correlazione è al di sopra di una certa soglia, poi andare avanti altrimenti rinunciare e andare avanti.

Si potrebbe avere un "controllare ogni 8 ° posizione" che si moltiplica per 8. Se questo è troppo basso, saltare e andare avanti. Se questo è abbastanza alto, quindi controllare tutti i valori per vedere se hai trovato i massimi.

Il problema è il tempo necessario per fare tutte queste moltiplica - (f[<unsigned int>(i+j)] * g[j]) In effetti, si sta riempiendo una grande matrice con tutti questi prodotti e raccogliendo la riga con la somma massima. Non si vuole calcolare "tutti" i prodotti. Appena sufficiente dei prodotti per essere sicuri di aver trovato l'importo massimo.

Il problema con la ricerca di massimi è che si deve sommare tutto per vedere se è più grande. Se si riesce a trasformare questo in un problema di minimizzazione, è più facile abbandonare i prodotti e le somme di calcolo una volta un risultato intermedio supera una soglia.

(credo che questo potrebbe funzionare. Ho have't provato.)

Se è stato utilizzato max(g)-g[j] di lavorare con numeri negativi, devi essere alla ricerca per i più piccoli, non il più grande. Si potrebbe calcolare la correlazione per la prima posizione. Tutto ciò che ha riassunto in un valore più grande poteva essere fermato immediatamente -. Niente più si moltiplica o aggiunge a quella offset, passare ad un altro

Answer 3

• è possibile estrarre gamma (size2) dal circuito esterno
• è possibile utilizzare sum () al posto di un ciclo per calcolare current_correlation
• è possibile memorizzare le correlazioni e ritardi in un elenco e quindi utilizzare max () per ottenere la più grande