-
21-12-2019 - |
题
我正在尝试检测我的回声 叽叽喳喳 在我的声音中 记录 在 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]))
我根本不明白,但似乎很容易实现。这就是说我失败了(假设我的 修正1 是正确的)。我感觉我完全误解了这一点?
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)?或者我生成数组来测试它们的任意方式?或者我应该进行共轭而不是反转它?也就是说,性能并不是真正的问题,因此除非有明显的问题,否则您不必费心指出优化。
解决方案
您的实施 xcorr1
确实符合互相关的标准信号处理定义。
相对于您在开头添加零的询问:添加 chirp.length-1
零将使结果的索引 0 对应于传输的开始。但请注意,相关输出的峰值出现 chirp.length-1
回声开始后的样本(线性调频脉冲必须与完整接收到的回声对齐)。使用峰值索引来获取回波延迟,然后您必须通过减去延迟或丢弃第一个延迟来调整相关器延迟 chirp.length-1
输出结果。请注意,附加的零对应于开头的许多额外输出,因此您最好不要首先在开头添加这些零。
为了 xcorr2
然而,有一些事情需要解决。首先,如果 recording
和 chirp
输入尚未被零填充至至少 chirp+recording 数据 您需要这样做的长度(出于性能原因,最好是 2 的幂长度)。如您所知,它们都需要填充到相同的长度。
其次,您没有考虑到中所示的乘法 发布参考答案, ,实际上对应于复杂的乘法(而 DoubleFFT_1D.realForward
API 使用双精度)。现在,如果您要使用线性调频脉冲的 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 的幂)。
其他提示
直接版本不需要首先要求零填充。您只需录制长度生成的长度生成播放器和啁啾的长度生成的曲线,并计算长度生成的结果。通过手工通过一个小的例子进行工作以进行步骤:
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输出是复杂的数字,因此您需要使用复杂乘法 。 (1 + 2J)*(3 + 4J)为-5 + 10J,而不是3 + 8J。我不知道复杂的数字是如何安排或处理的,但确保这是正确的。
或者是转换的目的;将实数移动到复杂域名?
否,傅立叶变换从时域变换到频域。时域数据可以是真实的或复杂的,并且频域数据可以是真实的或复杂的。在大多数情况下,您可以使用复杂的频谱具有真实数据。您需要在傅立叶变换上读取。
如果遥犯实际上是返回虚构/复杂的数字,那么复杂的情况如何?
真正的fft需要一个真正的输入,而复杂的fft需要一个复杂的输入。两个转换都会产生复杂的数字作为其输出。这就是DFT所做的。 DFT产生真实输出的唯一时间是输入数据是对称的(在这种情况下,您可以使用DCT来节省更多时间)。