2

我正在尝试在 Android 上的录音中检测我的啁啾回声，似乎互相关是找到两个信号的 FFT 相似位置的最合适方法，从那里我可以识别互相关数组中的峰值将对应于距离。

据我了解，我提出了以下互相关函数。这个对吗？我不确定是否在开头添加零并重新开始一些元素？

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 版本的两个数组的两端添加零。

在 SleuthEye 的回应后编辑：

您能否验证一下，因为我正在处理“实际”数据，所以我只需要通过进行真正的转换来完成一半的计算（实际部分）吗？

从您的代码中，看起来 REAL 转换返回的数组中的奇数元素是虚构的。这里发生了什么？

我如何从实数数组变为复数？或者这是转换的目的吗？将实数移入复数域？（但实数只是复数的一个子集，所以它们不是已经在这个域中了吗？）

如果 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）？或者我生成数组来测试它们的任意方式？或者我应该做共轭而不是逆转它？也就是说，性能并不是真正的问题，因此除非有明显的东西，否则您无需费心指出优化。

4

2 回答 2

2

您的实现xcorr1确实对应于互相关的标准信号处理定义。

相对于您关于在开头添加零的询问：添加chirp.length-1零将使结果的索引 0 对应于传输的开始。但是请注意，相关输出的峰值出现chirp.length-1在回波开始之后的样本（啁啾必须与完整接收到的回波对齐）。使用峰值指数来获得回波延迟，然后您必须通过减去延迟或丢弃第一个chirp.length-1输出结果来调整该相关器延迟。请注意，额外的零对应于开头的许多额外输出，您最好不要在开头添加这些零。

然而，有几xcorr2件事需要解决。首先，如果recording和chirp输入还没有被零填充到至少啁啾+记录数据长度，你需要这样做（出于性能原因，最好是长度的 2 次方）。如您所知，它们都需要填充到相同的长度。

其次，您没有考虑到发布的参考答案中指出的乘法实际上对应于复数乘法（而DoubleFFT_1D.realForwardAPI 使用双精度数）。现在，如果您要实现诸如与啁啾的 FFT 的复数乘法之类的东西，您还不如实际实现与啁啾 FFT 的复共轭的乘法（参考答案中指示的替代实现），无需反转时域值。

还考虑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 的幂）。

于 2014-07-31T16:37:09.170 回答
1

直接版本不需要首先填充零。您只需记录长度M和长度的啁啾N并计算长度的结果N+M-1。手动完成一个小例子来了解这些步骤：

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。我不知道您的复数是如何排列或处理的，但请确保这是正确的。

或者这是转换的目的吗？将实数移入复数域？

不，傅里叶变换从时域变换到频域。时域数据可以是实数或复数，频域数据可以是实数或复数。在大多数情况下，您拥有具有复杂光谱的真实数据。您需要阅读傅里叶变换。

如果 realForward 实际上返回的是虚数/复数，它与 complexForward 有何不同？

实数 FFT 采用实数输入，而复数 FFT 采用复数输入。两种转换都产生复数作为其输出。这就是 DFT 所做的。DFT 产生实际输出的唯一时间是输入数据是对称的（在这种情况下，您可以使用 DCT 来节省更多时间）。

于 2014-07-31T22:56:25.753 回答