Correlação cruzada para encontrar ecos de sonar

Question

Estou tentando detectar ecos do meu chilro no meu som gravação no Android e parece que a correlação cruzada é a maneira mais apropriada de descobrir onde as FFTs dos dois sinais são semelhantes e a partir daí posso identificar picos na matriz correlacionada cruzada que corresponderá às distâncias.

Pelo que entendi, criei a seguinte função de correlação cruzada.Isso está correto?Eu não tinha certeza se deveria adicionar zeros ao início e começar alguns 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;
}

Pergunta secundária:

De acordo com esse resposta, também pode ser calculado como

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

que eu não entendo nada, mas parece fácil de implementar.Dito isto, falhei (assumindo que meu xcorr1 está correto).Eu sinto que entendi isso completamente mal?

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

Supondo que ambos funcionem, qual função seria mais apropriada, visto que os arrays conterão alguns milhares de elementos?

EDITAR:A propósito, tentei adicionar zeros em ambas as extremidades de ambos os arrays para a versão FFT.

---

EDIT após a resposta do SleuthEye:

Você pode apenas verificar que, como estou lidando com dados “reais”, preciso fazer apenas metade dos cálculos (as partes reais) fazendo uma transformação real?

Pelo seu código, parece que os elementos ímpares na matriz retornada pela transformação REAL são imaginários.O que está acontecendo aqui?

Como vou de uma matriz de números reais para complexos?Ou é este o propósito de uma transformação;mover números reais para o domínio complexo?(mas os números reais são apenas um subconjunto dos números complexos e, portanto, já não estariam neste domínio?)

Se realForward está de fato retornando números imaginários/complexos, como ele difere de complexForward?E como interpreto os resultados?A magnitude do número complexo?

Peço desculpas pela minha falta de compreensão em relação às transformações, só até agora estudei as séries de Fourier.

Obrigado pelo código.Aqui está a 'minha' implementação de trabalho:

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

No entanto, até cerca de 5.000 elementos cada, xcorr1 parece ser mais rápido.Estou fazendo algo particularmente lento (talvez a constante atualização da memória - talvez eu deva converter para um ArrayList)?Ou a maneira arbitrária como gerei os arrays para testá-los?Ou devo fazer os conjugados em vez de reverter?Dito isto, o desempenho não é realmente um problema, então, a menos que haja algo óbvio, você não precisa se preocupar em apontar otimizações.

Answer 1

Sua implementação de xcorr1 corresponde à definição padrão de correlação cruzada de processamento de sinal.

Em relação à sua interrogação em relação à adição de zeros no início:adicionando chirp.length-1 zeros fariam com que o índice 0 do resultado correspondesse ao início da transmissão.Observe, entretanto, que o pico da saída da correlação ocorre chirp.length-1 amostras após o início dos ecos (o chirp deve estar alinhado com o eco recebido completo).Usando o índice de pico para obter atrasos de eco, você teria que ajustar o atraso do correlacionador subtraindo o atraso ou descartando o primeiro chirp.length-1 resultados de saída.Observando que os zeros adicionais correspondem a muitas saídas extras no início, provavelmente seria melhor não adicionar esses zeros no início.

Para xcorr2 no entanto, algumas coisas precisam ser abordadas.Primeiro, se o recording e chirp as entradas ainda não foram preenchidas com zeros para pelo menos chirp + gravação dados comprimento, você precisaria fazer isso (de preferência com uma potência de 2 comprimento por motivos de desempenho).Como você sabe, ambos precisariam ser preenchidos no mesmo comprimento.

Segundo, você não levou em conta que a multiplicação indicada no resposta de referência postada, correspondem de fato a multiplicações complexas (enquanto DoubleFFT_1D.realForward API usa duplas).Agora, se você for implementar algo como uma multiplicação complexa com a FFT do chirp, você também pode implementar a multiplicação com o conjugado complexo da FFT do chirp (a implementação alternativa indicada no resposta de referência), eliminando a necessidade de reverter os valores no domínio do tempo.

Também contabilizando DoubleFFT_1D.realForward ordem de embalagem para transformações de comprimento par, você obteria:

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

Observe que o result matriz seria do mesmo tamanho que paddedRecording e paddedChirp, mas apenas o primeiro recording.length+chirp.length-1 deve ser mantido.

Finalmente, em relação a qual função é mais apropriada para arrays de alguns milhares de elementos, a versão FFT xcorr2 provavelmente será muito mais rápido (desde que você restrinja os comprimentos da matriz a potências de 2).

Answer 2

A versão direta não requer preenchimento de zeros primeiro.Você apenas grava a duração M e chilrear de comprimento N e calcular o resultado do comprimento N+M-1.Trabalhe com um pequeno exemplo manualmente para entender as etapas:

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]

O método FFT é muito mais rápido se você tiver arrays longos.Neste caso, você deve preencher cada entrada com zero no tamanho M + N-1 para que ambas as matrizes de entrada tenham o mesmo tamanho antes tomando a FFT.

Além disso, a saída FFT contém números complexos, então você precisa usar multiplicação complexa.(1+2j)*(3+4j) é -5+10j, não 3+8j.Não sei como seus números complexos são organizados ou tratados, mas certifique-se de que isso esteja correto.

Ou é este o propósito de uma transformação;mover números reais para o domínio complexo?

Não, a transformada de Fourier transforma do domínio do tempo para o domínio da frequência.Os dados no domínio do tempo podem ser reais ou complexos, e os dados no domínio da frequência podem ser reais ou complexos.Na maioria dos casos você tem dados reais com um espectro complexo.Você precisa ler sobre a transformada de Fourier.

Se realForward está de fato retornando números imaginários/complexos, como ele difere de complexForward?

A verdadeira FFT leva um verdadeiro entrada, enquanto a FFT complexa assume um complexo entrada.Ambas as transformações produzem números complexos como saída.É isso que o DFT faz.A única vez que uma DFT produz uma saída real é se os dados de entrada forem simétricos (nesse caso, você pode usar a DCT para economizar ainda mais tempo).