android - 互相关以找到声纳回波

我试图在android上检测我的声音记录中chirp的回声，似乎互相关是找到两个信号的fft相似的地方的最合适的方法，从那里我可以识别出互相关阵列中与距离相对应的峰值。
根据我的理解，我提出了以下互相关函数。这是对的吗？我不确定是否要在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;
}

次要问题：
根据this答案，也可以计算如下

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）？或者我生成数组来测试它们的任意方式？或者我应该做共轭而不是反转它？也就是说，性能并不是一个真正的问题，所以除非有明显的东西，你不必费心指出优化。

最佳答案

您的xcorr1实现与互相关的标准信号处理定义相对应。
相对于在开始时加零的询问：加上chirp.length-1零将使结果的索引0对应于传输的开始。但是，请注意，相关输出的峰值在回波开始后出现chirp.length-1采样（啁啾必须与接收到的全部回波对齐）。使用峰值索引获得回波延迟，然后必须通过减去延迟或丢弃第一个chirp.length-1输出结果来调整相关器延迟。注意到额外的零对应于开始时的许多额外输出，您最好不要首先在开始时添加这些零。
然而，对于xcorr2来说，有一些事情需要解决。首先，如果recording和chirp输入尚未零填充到至少chirp+记录数据长度，则需要这样做（出于性能原因，最好是2倍长度的功率）。如你所知，它们都需要填充到相同的长度。
其次，您没有考虑到posted reference answer中表示的乘法实际上对应于复杂的乘法（而DoubleFFT_1D.realForwardapi使用双倍）。现在，如果您要实现与chirp的fft的复数乘法之类的东西，您也可以实际实现与chirp的fft的复数共轭（在reference answer中指示的替代实现）的乘法，从而不需要反转时域值。
同样考虑到偶数长度变换的打包顺序，您将得到：

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

注意DoubleFFT_1D.realForward数组的大小与result和paddedRecording相同，但只应保留第一个paddedChirp。
最后，相对于哪个函数最适合几千个元素的数组，fft版本recording.length+chirp.length-1可能要快得多（前提是将数组长度限制为2的幂）。