Kreuzkorrelation zum Auffinden von Sonar-Echos

Question

Ich versuche, Echos von mir zu erkennen zwitschern in meinem Klang Aufzeichnung auf Android und es scheint, dass Kreuzkorrelation die am besten geeignete Methode ist, um herauszufinden, wo die FFTs der beiden Signale ähnlich sind, und von dort aus kann ich Spitzen im kreuzkorrelierten Array identifizieren, die Entfernungen entsprechen.

Nach meinem Verständnis habe ich die folgende Kreuzkorrelationsfunktion entwickelt.Ist das richtig?Ich war mir nicht sicher, ob ich Nullen am Anfang hinzufügen und ein paar Elemente zurück beginnen sollte?

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

Sekundäre Frage:

Entsprechend Das Antwort, es kann auch wie berechnet werden

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

was ich überhaupt nicht verstehe, aber einfach genug zu implementieren scheint.Das heißt, ich habe versagt (vorausgesetzt, mein xcorr1 ist richtig).Ich habe das Gefühl, dass ich das völlig falsch verstanden habe?

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

Angenommen, ich habe beides zum Laufen gebracht. Welche Funktion wäre angesichts der Tatsache, dass die Arrays einige tausend Elemente enthalten, am besten geeignet?

BEARBEITEN:Übrigens habe ich versucht, für die FFT-Version an beiden Enden beider Arrays Nullen hinzuzufügen.

---

EDIT nach der Antwort von SleuthEye:

Können Sie einfach bestätigen, dass ich nur die Hälfte der Berechnungen (die realen Teile) durchführen muss, indem ich eine echte Transformation durchführe, da ich es mit „tatsächlichen“ Daten zu tun habe?

Aus Ihrem Code geht hervor, dass die ungeradzahligen Elemente im von der REAL-Transformation zurückgegebenen Array imaginär sind.Was ist denn hier los?

Wie komme ich von einem Array reeller Zahlen zu komplexen Zahlen?Oder ist dies der Zweck einer Transformation?reelle Zahlen in den komplexen Bereich verschieben?(Aber die reellen Zahlen sind nur eine Teilmenge der komplexen Zahlen und wären sie also nicht bereits in diesem Bereich?)

Wenn realForward tatsächlich imaginäre/komplexe Zahlen zurückgibt, worin besteht dann der Unterschied zu complexForward?Und wie interpretiere ich die Ergebnisse?Der Betrag der komplexen Zahl?

Ich entschuldige mich für mein mangelndes Verständnis in Bezug auf Transformationen, ich habe bisher nur Fourier-Reihen studiert.

Danke für den Code.Hier ist „meine“ funktionierende Implementierung:

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

Allerdings scheint xcorr1 bis jeweils etwa 5000 Elemente schneller zu sein.Mache ich etwas besonders Langsames (vielleicht das ständige „Neuerstellen“ des Speichers – vielleicht sollte ich in eine ArrayList umwandeln)?Oder die willkürliche Art und Weise, wie ich die Arrays generiert habe, um sie zu testen?Oder sollte ich die Konjugate machen, anstatt es umzukehren?Allerdings ist die Leistung kein wirkliches Problem, sodass Sie sich nicht die Mühe machen müssen, auf Optimierungen hinzuweisen, es sei denn, es liegt etwas Offensichtliches vor.

Answer 1

Ihre Umsetzung von xcorr1 entspricht der Standard-Signalverarbeitungsdefinition der Kreuzkorrelation.

Bezogen auf Ihre Frage bezüglich des Hinzufügens von Nullen am Anfang:hinzufügen chirp.length-1 Nullen würden dazu führen, dass der Index 0 des Ergebnisses dem Beginn der Übertragung entspricht.Beachten Sie jedoch, dass der Höhepunkt der Korrelationsausgabe auftritt chirp.length-1 Abtastwerte nach Beginn der Echos (der Chirp muss mit dem gesamten empfangenen Echo übereinstimmen).Wenn Sie den Spitzenindex verwenden, um Echoverzögerungen zu erhalten, müssten Sie diese Korrelatorverzögerung dann entweder durch Subtrahieren der Verzögerung oder durch Verwerfen der ersten Verzögerung anpassen chirp.length-1 Ergebnisse ausgeben.Da die zusätzlichen Nullen so vielen zusätzlichen Ausgaben am Anfang entsprechen, wäre es wahrscheinlich besser, diese Nullen überhaupt nicht am Anfang hinzuzufügen.

Für xcorr2 Es müssen jedoch einige Dinge angegangen werden.Erstens, wenn die recording Und chirp Eingaben sind nicht bereits mit Nullen aufgefüllt, um mindestens Chirp+Aufzeichnung zu ermöglichen Daten Länge, die Sie dazu benötigen würden (aus Leistungsgründen vorzugsweise auf eine Potenz von 2 Länge).Wie Sie wissen, müssten beide auf die gleiche Länge gepolstert werden.

Zweitens haben Sie nicht berücksichtigt, dass die in der angegebenen Multiplikation Gepostete Referenzantwort, entsprechen in der Tat komplexen Multiplikationen (während DoubleFFT_1D.realForward API verwendet Doubles).Wenn Sie nun beispielsweise eine komplexe Multiplikation mit der FFT des Chirps implementieren möchten, können Sie die Multiplikation auch tatsächlich mit der komplex konjugierten FFT des Chirps implementieren (die in der Abbildung angegebene alternative Implementierung). Referenzantwort), wodurch die Notwendigkeit entfällt, die Zeitbereichswerte umzukehren.

Auch buchhalterisch DoubleFFT_1D.realForward Packreihenfolge für Transformationen mit gerader Länge erhalten Sie:

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

Notiere dass der result Das Array hätte die gleiche Größe wie paddedRecording Und paddedChirp, aber nur der erste recording.length+chirp.length-1 sollte zurückgehalten werden.

Schließlich ist die FFT-Version relativ dazu, welche Funktion für Arrays mit einigen tausend Elementen am besten geeignet ist xcorr2 wird wahrscheinlich viel schneller sein (vorausgesetzt, Sie beschränken die Array-Längen auf Potenzen von 2).

Answer 2

Die Direktversion erfordert zuerst keine Null-Polsterung. Sie nehmen einfach die Aufzeichnung von Längengroß-oditicetagcode und Chirp aus dem Längengroß-odizidagcode und berechnen das Ergebnis der Längengroß-odizidagcode. Arbeiten Sie durch ein winziges Beispiel von Hand, um die Schritte:

generasacodicetagpre.

Die FFT-Methode ist viel schneller, wenn Sie lange Arrays haben. In diesem Fall müssen Sie jeden Eingang mit jedem Eingang M + N-1 null ankissen, so dass beide Eingabearrays die gleiche Größe sind, bevor die FFT nimmt.

Außerdem ist der FFT-Ausgang komplexe Zahlen, sodass Sie komplexe Multiplikation verwenden müssen . (1 + 2j) * (3 + 4J) beträgt -5 + 10J, nicht 3 + 8J. Ich weiß nicht, wie Ihre komplexen Zahlen arrangiert oder gehandhabt werden, aber sicherstellen, dass dies richtig ist.

oder ist dies der Zweck einer Transformation; Um reelle Zahlen in die komplexe Domäne zu verschieben?

Nein, die Fourier-Transformation wandelt sich aus dem Zeitbereich in den Frequenzbereich um. Die Zeitdomänendaten können entweder echt oder komplex sein, und die Frequenzdomänendaten können entweder echt oder komplex sein. In den meisten Fällen haben Sie echte Daten mit einem komplexen Spektrum. Sie müssen auf der Fourier-Transformation überlesen.

Wenn revorward tatsächlich imaginäre / komplexe Zahlen zurückkehrt, wie unterscheidet es sich um komplexe Förderung?

Die echte FFT nimmt einen echten input , während die komplexe FFT einen komplexen input nimmt. Beide Transformationen erzeugen komplexe Zahlen als ihre Leistung. Das macht die DFT. Die einzige Zeit, in der ein DFT eine echte Ausgabe erzeugt, ist, wenn die Eingabedaten symmetrisch sind (In diesem Fall können Sie das DCT verwenden, um noch mehr Zeit zu sparen).