Перекрестная корреляция для поиска эхо-сигналов гидролокатора

Question

Я пытаюсь уловить отголоски моего чириканье в моем звуке запись на Android, и, похоже, кросс-корреляция является наиболее подходящим способом определения того, где БПФ двух сигналов схожи, и оттуда я могу идентифицировать пики в кросс-коррелированном массиве, которые будут соответствовать расстояниям.

Насколько я понимаю, я получил следующую функцию взаимной корреляции.Правильно ли это?Я не был уверен, стоит ли добавлять нули в начало as и начинать несколько элементов обратно?

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

Второстепенный вопрос:

Согласно это ответ, он также может быть рассчитан следующим образом

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

который я вообще не понимаю, но кажется достаточно простым в реализации.Тем не менее, я потерпел неудачу (предполагая, что мой xcorr1 является правильным).У меня такое чувство, что я совершенно неправильно это понял?

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

Предполагая, что у меня работают оба варианта, какая функция была бы наиболее подходящей, учитывая, что массивы будут содержать несколько тысяч элементов?

Редактировать:Кстати, я попытался добавить нули на оба конца обоих массивов для версии FFT.

---

РЕДАКТИРОВАТЬ после ответа сыщика:

Можете ли вы просто подтвердить, что, поскольку я имею дело с "фактическими" данными, мне нужно выполнить только половину вычислений (реальные части), выполнив реальное преобразование?

Из вашего кода это выглядит так, как будто элементы с нечетными номерами в массиве, возвращаемом РЕАЛЬНЫМ преобразованием, являются воображаемыми.Что здесь происходит?

Как мне перейти от массива действительных чисел к комплексному?Или это и есть цель преобразования;чтобы переместить действительные числа в сложную область?(но действительные числа - это всего лишь подмножество комплексных чисел, и поэтому разве они уже не были бы в этой области?)

Если realForward на самом деле возвращает мнимые / комплексные числа, чем это отличается от complexForward?И как мне интерпретировать полученные результаты?Величина комплексного числа?

Я приношу извинения за мое непонимание относительно преобразований, я пока только изучал ряды Фурье.

Спасибо за код.Вот "моя" рабочая реализация:

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

Однако, примерно до 5000 элементов в каждом, xcorr1, кажется, работает быстрее.Делаю ли я что-нибудь особенно медленно (возможно, постоянное "обновление" памяти - может быть, мне следует привести к ArrayList)?Или произвольный способ, которым я сгенерировал массивы для их тестирования?Или я должен делать конъюгаты вместо того, чтобы менять их местами?Тем не менее, производительность на самом деле не является проблемой, поэтому, если нет чего-то очевидного, вам не нужно утруждать себя указанием на оптимизацию.

Answer 1

Ваша реализация xcorr1 соответствует стандартному определению взаимной корреляции при обработке сигналов.

Относительно вашего запроса относительно добавления нулей в начале:добавление chirp.length-1 из-за нулей индекс 0 результата будет соответствовать началу передачи.Однако обратите внимание, что происходит пик выходного сигнала корреляции chirp.length-1 выборки после начала эхо-сигналов (чириканье должно быть выровнено с полным полученным эхо-сигналом).Используя пиковый индекс для получения задержек эхо-сигнала, вам затем придется скорректировать задержку этого коррелятора либо путем вычитания задержки, либо путем отбрасывания первой chirp.length-1 выводите результаты.Отмечая, что дополнительные нули соответствуют такому количеству дополнительных выходных данных в начале, вам, вероятно, было бы лучше вообще не добавлять эти нули в начале.

Для xcorr2 однако необходимо решить несколько вопросов.Во-первых, если recording и chirp входные данные еще не обнулены, по крайней мере, для записи chirp + данные длина, которую вам нужно было бы сделать таким образом (предпочтительно в степени 2 длины по соображениям производительности).Как вы знаете, они оба должны быть дополнены до одинаковой длины.

Во-вторых, вы не приняли во внимание, что умножение, указанное в опубликованный ссылочный ответ, фактически соответствуют комплексным умножениям (тогда как DoubleFFT_1D.realForward API использует дубли).Теперь, если вы собираетесь реализовать что-то вроде комплексного умножения с помощью FFT chirp, вы могли бы также фактически реализовать умножение с комплексно сопряженным FFT chirp (альтернативная реализация, указанная в справочный ответ), устраняя необходимость изменять значения во временной области.

Также учитывается DoubleFFT_1D.realForward при заказе упаковки для преобразования четной длины вы получите:

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

Обратите внимание , что result массив был бы того же размера, что и paddedRecording и paddedChirp, но только первый recording.length+chirp.length-1 должно быть сохранено.

Наконец, относительно того, какая функция является наиболее подходящей для массивов из нескольких тысяч элементов, версия FFT xcorr2 вероятно, это будет намного быстрее (при условии, что вы ограничите длину массива степенями 2).

Answer 2

Прямая версия не требует от заполнения нулевой прокладки. Вы просто принимаете запись о длине Genacodicetacodcode и Chirp длины генеракодицетагкода и рассчитывают результат длины генеракодицетагкода. Работайте через крошечный пример вручную, чтобы разделить шаги:

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]

.

Метод FFT намного быстрее, если у вас есть длинные массивы. В этом случае у вас есть нулевой панель каждый вход на размер M + N-1, чтобы оба входных массива были одинаковыми размерами до , принимая FFT.

Кроме того, вывод FFT - это сложные числа, поэтому вам нужно использовать Комплексное умножение Отказ (1 + 2J) * (3 + 4J) - -5 + 10J, а не 3 + 8J. Я не знаю, как ваши сложные числа расположены или обрабатываются, но убедитесь, что это правильно.

или это цель преобразования; Переместить реальные номера в комплексный домен?

Нет, преобразование Фурье трансформируется из временного домена к частотному домену. Данные временного домена могут быть реальными или сложными, и данные частоты домена могут быть либо реальными, либо сложными. В большинстве случаев у вас есть реальные данные со сложным спектром. Вы должны прочитать в преобразовании Фурье.

Если реализуется на самом деле возвращение воображаемых / комплексных чисел, как он отличается от ComperyForward?

Настоящий FFT принимает реальный <прочный> вход , в то время как комплекс FFT принимает комплекс вход . Оба преобразования производят комплексные числа в качестве их выхода. Это то, что делает ДФТ. Единственный раз, когда DFT производит реальный выход, заключается в том, что входные данные являются симметричными (в этом случае вы можете использовать DCT для сохранения еще больше времени).