Correlación cruzada para encontrar ecos de sonar

Question

Estoy tratando de detectar ecos de mi chirrido en mi sonido grabación en Android y parece que la correlación cruzada es la forma más adecuada de encontrar dónde son similares las FFT de las dos señales y a partir de ahí puedo identificar picos en la matriz de correlación cruzada que corresponderán a las distancias.

Según tengo entendido, se me ocurrió la siguiente función de correlación cruzada.¿Es esto correcto?¿No estaba seguro de si agregar ceros al principio y comenzar algunos elementos atrás?

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;
}

Pregunta secundaria:

De acuerdo a este respuesta, también se puede calcular como

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

Lo cual no entiendo en absoluto pero parece bastante fácil de implementar.Dicho esto, he fallado (asumiendo mi xcorr1 es correcto).¿Siento que he entendido mal esto por completo?

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;
}

Suponiendo que ambas funcionen, ¿qué función sería la más apropiada dado que las matrices contendrán unos pocos miles de elementos?

EDITAR:Por cierto, intenté agregar ceros a ambos extremos de ambas matrices para la versión FFT.

---

EDITAR después de la respuesta de SleuthEye:

¿Puedes simplemente verificar que, debido a que estoy tratando con datos "reales", solo necesito hacer la mitad de los cálculos (las partes reales) mediante una transformación real?

Según su código, parece que los elementos impares en la matriz devuelta por la transformación REAL son imaginarios.¿Que está pasando aqui?

¿Cómo paso de una serie de números reales a complejos?¿O es este el propósito de una transformación?¿Mover números reales al dominio complejo?(pero los números reales son solo un subconjunto de los números complejos y, por lo tanto, ¿no estarían ya en este dominio?)

Si realForward de hecho devuelve números imaginarios/complejos, ¿en qué se diferencia de complexForward?¿Y cómo interpreto los resultados?¿La magnitud del número complejo?

Pido disculpas por mi falta de comprensión con respecto a las transformadas, hasta ahora solo he estudiado series de Fourier.

Gracias por el código.Aquí está "mi" implementación de trabajo:

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;
}

Sin embargo, hasta aproximadamente 5000 elementos cada uno, xcorr1 parece ser más rápido.¿Estoy haciendo algo particularmente lento (quizás la constante "nueva" memoria, tal vez debería convertir una ArrayList)?¿O la forma arbitraria en la que generé los arreglos para probarlos?¿O debería hacer los conjugados en lugar de invertirlo?Dicho esto, el rendimiento no es realmente un problema, por lo que, a menos que haya algo obvio, no es necesario molestarse en señalar optimizaciones.

Answer 1

Su implementación de xcorr1 corresponde a la definición estándar de correlación cruzada del procesamiento de señales.

En relación con su interrogatorio con respecto a agregar ceros al principio:agregando chirp.length-1 los ceros harían que el índice 0 del resultado corresponda al inicio de la transmisión.Sin embargo, tenga en cuenta que el pico de la salida de correlación se produce chirp.length-1 muestras después del inicio de los ecos (el chirrido debe estar alineado con el eco recibido completo).Usando el índice de pico para obtener retrasos del eco, entonces tendría que ajustar ese retraso del correlador, ya sea restando el retraso o descartando el primer chirp.length-1 resultados de salida.Teniendo en cuenta que los ceros adicionales corresponden a tantas salidas adicionales al principio, probablemente sería mejor no agregar esos ceros al principio en primer lugar.

Para xcorr2 sin embargo, es necesario abordar algunas cosas.Primero, si el recording y chirp las entradas aún no están rellenadas con ceros para al menos chirrido + grabación datos longitud que necesitaría para hacerlo (preferiblemente a una potencia de 2 longitud por razones de rendimiento).Como ya sabes, ambos deberían estar acolchados con la misma longitud.

Segundo, no tomaste en cuenta que la multiplicación indicada en el respuesta de referencia publicada, corresponden de hecho a multiplicaciones complejas (mientras que DoubleFFT_1D.realForward API utiliza dobles).Ahora bien, si va a implementar algo como una multiplicación compleja con la FFT del chirrido, también podría implementar la multiplicación con el conjugado complejo de la FFT del chirrido (la implementación alternativa indicada en el respuesta de referencia), eliminando la necesidad de invertir los valores en el dominio del tiempo.

También contabilizando DoubleFFT_1D.realForward orden de embalaje para transformaciones de longitud uniforme, obtendría:

// [...]
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);
// [...]

Tenga en cuenta que el result la matriz sería del mismo tamaño que paddedRecording y paddedChirp, pero sólo el primero recording.length+chirp.length-1 debería ser guardado.

Finalmente, en relación con qué función es la más apropiada para matrices de unos pocos miles de elementos, la versión FFT xcorr2 Es probable que sea mucho más rápido (siempre que restrinja la longitud de la matriz a potencias de 2).

Answer 2

La versión directa no requiere relleno con ceros primero.Solo tomas la grabación de la duración. M y chirrido de longitud N y calcular el resultado de la longitud N+M-1.Trabaje a mano con un pequeño ejemplo para asimilar los pasos:

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]

El método FFT es mucho más rápido si tiene matrices largas.En este caso, debe rellenar con ceros cada entrada al tamaño M+N-1 para que ambas matrices de entrada tengan el mismo tamaño. antes tomando la FFT.

Además, la salida FFT son números complejos, por lo que es necesario utilizar multiplicación compleja.(1+2j)*(3+4j) es -5+10j, no 3+8j.No sé cómo se organizan o manejan sus números complejos, pero asegúrese de que sea correcto.

¿O es este el propósito de una transformación?¿Mover números reales al dominio complejo?

No, la transformada de Fourier pasa del dominio del tiempo al dominio de la frecuencia.Los datos en el dominio del tiempo pueden ser reales o complejos, y los datos en el dominio de la frecuencia pueden ser reales o complejos.En la mayoría de los casos se tienen datos reales con un espectro complejo.Necesita leer más sobre la transformada de Fourier.

Si realForward de hecho devuelve números imaginarios/complejos, ¿en qué se diferencia de complexForward?

La FFT real toma un valor real aporte, mientras que la FFT compleja toma una compleja aporte.Ambas transformaciones producen números complejos como resultado.Eso es lo que hace el DFT.La única vez que una DFT produce una salida real es si los datos de entrada son simétricos (en cuyo caso puede usar la DCT para ahorrar aún más tiempo).