Corrélation croisée pour trouver les échos du sonar

Question

J'essaie de détecter les échos de mon gazouiller dans mon son enregistrement sur Android et il semble que la corrélation croisée soit le moyen le plus approprié de trouver où les FFT des deux signaux sont similaires et à partir de là, je peux identifier les pics dans le tableau à corrélation croisée qui correspondront aux distances.

D'après ma compréhension, j'ai proposé la fonction de corrélation croisée suivante.Est-ce correct?Je ne savais pas s'il fallait ajouter des zéros au début et recommencer quelques éléments ?

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

Question secondaire :

Selon ce réponse, cela peut aussi être calculé comme

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

ce que je ne comprends pas du tout mais qui semble assez simple à mettre en œuvre.Cela dit, j'ai échoué (en supposant que mon xcorr1 est correct).J'ai l'impression d'avoir complètement mal compris ?

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

En supposant que les deux fonctionnent, quelle fonction serait la plus appropriée étant donné que les tableaux contiendront quelques milliers d'éléments ?

MODIFIER:Au fait, j'ai essayé d'ajouter des zéros aux deux extrémités des deux tableaux pour la version FFT.

---

EDIT après la réponse de SleuthEye :

Pouvez-vous simplement vérifier que, parce que je traite des données « réelles », je n'ai besoin de faire que la moitié des calculs (les parties réelles) en effectuant une véritable transformation ?

D'après votre code, il semble que les éléments impairs du tableau renvoyés par la transformation REAL soient imaginaires.Que se passe t-il ici?

Comment puis-je passer d’un tableau de nombres réels à des nombres complexes ?Ou est-ce le but d’une transformation ?déplacer les nombres réels dans le domaine complexe ?(mais les nombres réels ne sont qu'un sous-ensemble des nombres complexes et ne seraient-ils donc pas déjà dans ce domaine ?)

Si realForward renvoie en fait des nombres imaginaires/complexes, en quoi diffère-t-il de complexForward ?Et comment interpréter les résultats ?La grandeur du nombre complexe ?

Je m'excuse pour mon manque de compréhension en ce qui concerne les transformées, je n'ai jusqu'à présent étudié que les séries de Fourier.

Merci pour le code.Voici « ma » mise en œuvre fonctionnelle :

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

Cependant, jusqu'à environ 5 000 éléments chacun, xcorr1 semble être plus rapide.Est-ce que je fais quelque chose de particulièrement lent (peut-être le « renouvellement » constant de la mémoire – peut-être devrais-je effectuer un cast vers une ArrayList) ?Ou la manière arbitraire dont j'ai généré les tableaux pour les tester ?Ou dois-je faire les conjugués au lieu de l'inverser ?Cela dit, les performances ne sont pas vraiment un problème, donc à moins qu'il n'y ait quelque chose d'évident, vous n'avez pas besoin de signaler les optimisations.

Answer 1

Votre mise en œuvre de xcorr1 correspond à la définition standard de traitement du signal de la corrélation croisée.

Par rapport à votre interrogation concernant l'ajout de zéros au début :ajouter chirp.length-1 des zéros feraient correspondre l'index 0 du résultat au début de la transmission.Notez cependant que le pic de la sortie de corrélation se produit chirp.length-1 échantillons après le début des échos (le gazouillis doit être aligné avec l’écho complet reçu).En utilisant l'indice de crête pour obtenir les délais d'écho, vous devrez alors ajuster ce délai de corrélateur soit en soustrayant le délai, soit en supprimant le premier. chirp.length-1 résultats de sortie.En notant que les zéros supplémentaires correspondent à autant de sorties supplémentaires au début, vous feriez probablement mieux de ne pas ajouter ces zéros au début.

Pour xcorr2 cependant, quelques points doivent être réglés.Premièrement, si le recording et chirp les entrées ne sont pas déjà complétées par des zéros pour au moins chirp + enregistrement données longueur dont vous auriez besoin pour le faire (de préférence à une puissance de 2 pour des raisons de performances).Comme vous le savez, ils devront tous deux être rembourrés à la même longueur.

Deuxièmement, vous n'avez pas tenu compte du fait que la multiplication indiquée dans le réponse de référence publiée, correspondent en fait à des multiplications complexes (alors que DoubleFFT_1D.realForward L'API utilise des doubles).Maintenant, si vous envisagez d'implémenter quelque chose comme une multiplication complexe avec la FFT du chirp, vous pourriez tout aussi bien implémenter la multiplication avec le conjugué complexe de la FFT du chirp (l'implémentation alternative indiquée dans le réponse de référence), supprimant la nécessité d'inverser les valeurs du domaine temporel.

Comptabilisant également DoubleFFT_1D.realForward ordre d'emballage pour des transformations de longueur paire, vous obtiendriez :

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

Notez que le result le tableau serait de la même taille que paddedRecording et paddedChirp, mais seulement le premier recording.length+chirp.length-1 Devrait être gardé.

Enfin, relativement à la fonction la plus appropriée pour des tableaux de quelques milliers d'éléments, la version FFT xcorr2 sera probablement beaucoup plus rapide (à condition que vous limitiez la longueur des tableaux à des puissances de 2).

Answer 2

La version directe ne nécessite pas de remplissage nul au préalable.Vous prenez juste un enregistrement de durée M et gazouillis de longueur N et calculer le résultat de la longueur N+M-1.Travaillez à la main sur un petit exemple pour comprendre les étapes :

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]

La méthode FFT est beaucoup plus rapide si vous disposez de longs tableaux.Dans ce cas, vous devez mettre à zéro chaque entrée à la taille M+N-1 afin que les deux tableaux d'entrée aient la même taille. avant prendre la FFT.

De plus, la sortie FFT est constituée de nombres complexes, vous devez donc utiliser multiplication complexe.(1+2j)*(3+4j) vaut -5+10j, pas 3+8j.Je ne sais pas comment vos nombres complexes sont organisés ou traités, mais assurez-vous que c'est correct.

Ou est-ce le but d’une transformation ?déplacer les nombres réels dans le domaine complexe ?

Non, la transformée de Fourier passe du domaine temporel au domaine fréquentiel.Les données du domaine temporel peuvent être réelles ou complexes, et les données du domaine fréquentiel peuvent être réelles ou complexes.Dans la plupart des cas, vous disposez de données réelles avec un spectre complexe.Vous devez vous renseigner sur la transformée de Fourier.

Si realForward renvoie en fait des nombres imaginaires/complexes, en quoi diffère-t-il de complexForward ?

La vraie FFT prend une vraie saisir, alors que la FFT complexe prend un complexe saisir.Les deux transformations produisent des nombres complexes en sortie.C'est ce que fait le DFT.La seule fois où un DFT produit une sortie réelle, c'est si les données d'entrée sont symétriques (auquel cas vous pouvez utiliser le DCT pour gagner encore plus de temps).