Correlazione incrociata per trovare echi sonar

Question

Sto cercando di rilevare gli echi dei miei chirp nel mio suono registrazione su Android e sembra che il cross correlazione sia il modo più appropriato per trovare dove le FFT del Due segnali sono simili e da lì posso identificare picchi nell'array croce che corrisponde a distanze.

Dalla mia comprensione, ho trovato la seguente funzione di correlazione incrociata. È corretto? Non ero sicuro se aggiungere Zeros all'inizio e iniziare alcuni elementi indietro?

public double[] xcorr1(double[] recording, double[] chirp) {        
    double[] recordingZeroPadded = new double[recording.length + chirp.length];

    for (int i = recording.length; i < recording.length + chirp.length; ++i)
            recordingZeroPadded[i] = 0;

    for (int i = 0; i < recording.length; ++i)
            recordingZeroPadded[i] = recording[i];

    double[] result = new double[recording.length + chirp.length - 1];

    for (int offset = 0; offset < recordingZeroPadded.length - chirp.length; ++offset)
        for (int i = 0; i < chirp.length; ++i)
            result[offset] += chirp[i] * recordingZeroPadded[offset + i];
    return result;
}

.

Domanda secondaria:

Secondo Questa risposta , può anche essere calcolata come

corr(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))

.

Che non capisco affatto ma sembra abbastanza facile da implementare. Detto quello che ho fallito (supponendo che i miei xcorr1 sia corretto). Mi sento come se avessi completamente frainteso questo?

public double[] xcorr2(double[] recording, double[] chirp) {
    // assume same length arguments for now
    DoubleFFT_1D fft = new DoubleFFT_1D(recording.length);
    fft.realForward(recording);
    reverse(chirp);
    fft.realForward(chirp);
    double[] result = new double[recording.length];

    for (int i = 0; i < result.length; ++i)
        result [i] = recording[i] * chirp[i];

    fft.realInverse(result, true);
    return result;
}

.

Supponendo che abbia avuto entrambi i lavori, quale funzione sarebbe più appropriata dato che gli array contengono alcuni migliaia di elementi?

Modifica: BTW, ho provato ad aggiungere zeri ad entrambe le estremità di entrambi gli array per la versione FFT.

---

.

Modifica dopo la risposta di Sleutheye:

Puoi semplicemente verificarlo, perché ho a che fare con "effettivi" dati, ho bisogno solo di fare la metà dei calcoli (le parti reali) facendo una trasformazione reale?

Dal tuo codice, sembra che gli elementi numerati dispari nell'array restituiti dalla trasformazione reale siano immaginari. Cosa sta succedendo qui?

Come vado da una serie di numeri reali a complessi? O è lo scopo di una trasformazione; Per spostare i numeri reali nel dominio complesso? (Ma i numeri reali sono solo un sottoinsieme dei numeri complessi e quindi non sarebbero già in questo dominio?)

Se in modo reale è in realtà il restituzione dei numeri immaginari / complessi, come differisce per complelire il complesso? E come interpretano i risultati? La grandezza del numero complesso?

Mi scuso per la mia mancanza di comprensione per quanto riguarda le trasformazioni, ho solo così tanto studiato serie Fourier.

Grazie per il codice. Ecco 'La mia "implementazione funzionante:

public double[] xcorr2(double[] recording, double[] chirp) {
    // pad to power of 2 for optimisation
    int y = 1;
    while (Math.pow(2,y) < recording.length + chirp.length)
        ++y;
    int paddedLength = (int)Math.pow(2,y);

    double[] paddedRecording = new double[paddedLength];
    double[] paddedChirp = new double[paddedLength];

    for (int i = 0; i < recording.length; ++i)
            paddedRecording[i] = recording[i];

    for (int i = recording.length; i < paddedLength; ++i)
            paddedRecording[i] = 0;

    for (int i = 0; i < chirp.length; ++i)
            paddedChirp[i] = chirp[i];

    for (int i = chirp.length; i < paddedLength; ++i)
            paddedChirp[i] = 0;

    reverse(chirp);
    DoubleFFT_1D fft = new DoubleFFT_1D(paddedLength);
    fft.realForward(paddedRecording);
    fft.realForward(paddedChirp);
    double[] result = new double[paddedLength];

    result[0] = paddedRecording[0] * paddedChirp[0]; // value at f=0Hz is real-valued
    result[1] = paddedRecording[1] * paddedChirp[1]; // value at f=fs/2 is real-valued and packed at index 1
    for (int i = 1; i < result.length / 2; ++i) {
        double a = paddedRecording[2*i];
        double b = paddedRecording[2*i + 1];
        double c = paddedChirp[2*i];
        double d = paddedChirp[2*i + 1];

        // (a+b*j)*(c-d*j) = (a*c+b*d) + (b*c-a*d)*j
        result[2*i]     = a*c + b*d;
        result[2*i + 1] = b*c - a*d;
    }

    fft.realInverse(result, true);

    // discard trailing zeros
    double[] result2 = new double[recording.length + chirp.length - 1];
    for (int i = 0; i < result2.length; ++i)
        result2[i] = result[i];

    return result2;
}

.

Tuttavia, fino a circa 5000 elementi ciascuno, XcorRR1 sembra essere più veloce. Sto facendo qualcosa di particolarmente lento (forse il nuovo "nuovo di memoria costante - forse dovrei lanciare ad un arrayList)? O il modo arbitrario in cui ho generato gli array per testarli? O dovrei fare i coniugati invece di reversirlo? Detto questo, le prestazioni non sono realmente un problema, quindi a meno che non ci sia qualcosa di ovvio che non devi preoccuparti di sottolineare gli ottimizzazione.

Answer 1

La tua implementazione di xcorr1 corrisponde alla definizione standard per l'elaborazione del segnale di correlazione cross-correlazione.

relativo al tuo interrogatorio rispetto all'aggiunta di zeri all'inizio: l'aggiunta di chirp.length-1 Zeros renderebbe l'indice 0 del risultato corrisponde all'inizio della trasmissione. Nota tuttavia che il picco dell'uscita di correlazione si verifica dei campioni chirp.length-1 dopo l'inizio degli echi (il chirp deve essere allineato con l'eco completo ricevuto). Utilizzando l'indice di picco per ottenere ritardi di eco, è necessario regolare il ritardo del correlatore sottraendo il ritardo o scartando i primi risultati di uscita chirp.length-1. Notando che gli zeri aggiuntivi corrispondono a quelle molte uscite aggiuntive all'inizio, probabilmente starai meglio di non aggiungere quegli zeri all'inizio in primo luogo.

Per xcorr2 Tuttavia, alcune cose devono essere affrontate. Innanzitutto, se gli ingressi recording e chirp non sono già a zero imbottiti per almeno chirp + registrazione dati è necessario farlo (preferibilmente a una potenza di 2 lunghezza per motivi di prestazione). Come sapete, entrambi dovrebbero essere imbottiti con la stessa lunghezza.

Secondo, non hai tenuto conto del fatto che la moltiplicazione indicata in risposta di riferimento pubblicato , corrisponde Fatto a moltiplicazioni complesse (mentre gli API DoubleFFT_1D.realForward utilizza doppi). Ora se hai intenzione di implementare qualcosa come una complessa moltiplicazione con il FFT del chirp, potrebbe anche implementare la moltiplicazione con il complesso coniugato della FFT del chirp (l'implementazione alternativa indicata in >> Risposta ), rimuovendo la necessità di invertire i valori del dominio del tempo.

Anche la contabilità per l'ordine di imballaggio DoubleFFT_1D.realForward per anche le trasformazioni di lunghezza, otterresti:

// [...]
fft.realForward(paddedRecording);
fft.realForward(paddedChirp);

result[0] = paddedRecording[0]*paddedChirp[0]; // value at f=0Hz is real-valued
result[1] = paddedRecording[1]*paddedChirp[1]; // value at f=fs/2 is real-valued and packed at index 1
for (int i = 1; i < result.length/2; ++i) {
    double a = paddedRecording[2*i];
    double b = paddedRecording[2*i+1];
    double c = paddedChirp[2*i];
    double d = paddedChirp[2*i+1];

    // (a+b*j)*(c-d*j) = (a*c+b*d) + (b*c-a*d)*j
    result[2*i]   = a*c + b*d;
    result[2*i+1] = b*c - a*d;
}
fft.realInverse(result, true);
// [...]

.

Notare che l'array di result avrebbe delle stesse dimensioni di paddedRecording e paddedChirp, ma deve essere mantenuta solo il primo recording.length+chirp.length-1.

Infine, rispetto a quale funzione è la funzione più appropriata per gli array di alcuni migliaia di elementi, è probabile che la versione FFT sia probabilmente molto più veloce (a condizione di limitare le lunghezze di array a poteri di 2).

Answer 2

La versione diretta non richiede prima il padding zero. Basta prendere la registrazione di lunghezza M e chirp di lunghezza N e calcola il risultato della lunghezza N+M-1. Lavorare attraverso un piccolo esempio a mano per aggrottare i passaggi:

recording = [1, 2, 3]
chirp = [4, 5]

  1 2 3
4 5

  1 2 3
  4 5

  1 2 3
    4 5

  1 2 3
      4 5


result = [1*5, 1*4 + 2*5, 2*4 + 3*5, 3*4] = [5, 14, 23, 4]

.

Il metodo FFT è molto più veloce se si dispongono di matrici lunghe. In questo caso è necessario eseguire il blocco a zero ciascun ingresso alla dimensione M + N-1 in modo che entrambi gli array di ingresso siano la stessa dimensione prima prendendo la FFT.

Anche l'uscita FFT è numeri complessi, quindi è necessario utilizzare moltiplicazione complessa . (1 + 2J) * (3 + 4J) è -5 + 10J, non 3 + 8J. Non so come i tuoi numeri complessi sono disposti o gestiti, ma assicurati che sia giusto.

.

o è lo scopo di una trasformazione; per spostare numeri reali nel dominio complesso?

No, la trasformazione di Fourier si trasforma dal dominio del tempo al dominio della frequenza. I dati del dominio del tempo possono essere reali o complessi, e i dati del dominio della frequenza possono essere reali o complessi. Nella maggior parte dei casi hai dati reali con uno spettro complesso. Devi leggere la trasformazione di Fourier.

.

Se in modo reale è in realtà il restituzione dei numeri immaginari / complessi, come differisce per complelire il complesso?

Il vero FFT prende un vero ingresso , mentre il complesso FFT prende un complesso Input . Entrambe le trasformazioni producono numeri complessi come loro output. Questo è ciò che fa il DFT. L'unica volta in cui un DFT produce un'uscita reale è se i dati di input sono simmetrici (nel qual caso è possibile utilizzare il DCT per risparmiare ancora più tempo).