寻找两个相似波形之间的时间偏移

Question

我需要比较两个时间和电压的波形图。由于这些波形的来源比较特殊，其中一个可能是另一个的时间偏移版本。

我该如何判断是否存在时间偏移？如果有的话，偏移量是多少呢？

我正在用Python来做这件事，并希望使用numpy和scipy这两个库。

Answer 1

这个函数对于真实值信号来说可能更有效。它使用了rfft，并且把输入数据填充为2的幂，这样可以确保进行线性（也就是非循环的）相关性计算：

def rfft_xcorr(x, y):
    M = len(x) + len(y) - 1
    N = 2 ** int(np.ceil(np.log2(M)))
    X = np.fft.rfft(x, N)
    Y = np.fft.rfft(y, N)
    cxy = np.fft.irfft(X * np.conj(Y))
    cxy = np.hstack((cxy[:len(x)], cxy[N-len(y)+1:]))
    return cxy

返回的值长度是 M = len(x) + len(y) - 1（通过 hstack 合并来去掉多余的零，因为我们要把长度凑成2的幂）。非负的延迟值是 cxy[0], cxy[1], ..., cxy[len(x)-1]，而负的延迟值是 cxy[-1], cxy[-2], ..., cxy[-len(y)+1]。

为了匹配一个参考信号，我会计算 rfft_xcorr(x, ref) 并寻找峰值。例如：

def match(x, ref):
    cxy = rfft_xcorr(x, ref)
    index = np.argmax(cxy)
    if index < len(x):
        return index
    else: # negative lag
        return index - len(cxy)   

In [1]: ref = np.array([1,2,3,4,5])
In [2]: x = np.hstack(([2,-3,9], 1.5 * ref, [0,3,8]))
In [3]: match(x, ref)
Out[3]: 3
In [4]: x = np.hstack((1.5 * ref, [0,3,8], [2,-3,-9]))
In [5]: match(x, ref)
Out[5]: 0
In [6]: x = np.hstack((1.5 * ref[1:], [0,3,8], [2,-3,-9,1]))
In [7]: match(x, ref)
Out[7]: -1

这并不是一个很稳妥的信号匹配方法，但它快速且简单。

Answer 2

如果一个信号相对于另一个信号有时间上的偏移，你会在它们的相关性中看到一个峰值。因为计算相关性需要消耗很多资源，所以用快速傅里叶变换（FFT）来处理会更好。因此，像下面这样的方法应该可以奏效：

af = scipy.fft(a)
bf = scipy.fft(b)
c = scipy.ifft(af * scipy.conj(bf))

time_shift = argmax(abs(c))

Answer 3

scipy提供了一个相关性函数，这个函数在处理小数据时效果很好。如果你想要非循环的相关性，也就是说信号不会绕回去，这个函数也能满足你的需求。需要注意的是，当你使用mode='full'时，signal.correlation返回的数组大小是两个信号大小的总和减去1（也就是len(a) + len(b) - 1），所以从argmax得到的值会比你预期的少一个信号的大小减去1（例如20）。

from scipy import signal, fftpack
import numpy
a = numpy.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0])
b = numpy.array([0, 0, 0, 0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0])
numpy.argmax(signal.correlate(a,b)) -> 16
numpy.argmax(signal.correlate(b,a)) -> 24

这两个不同的值对应于你是在看a的位移还是b的位移。

如果你想要循环相关性，并且信号比较大，可以使用卷积/傅里叶变换定理，但要注意，相关性和卷积非常相似，但并不完全相同。

A = fftpack.fft(a)
B = fftpack.fft(b)
Ar = -A.conjugate()
Br = -B.conjugate()
numpy.argmax(numpy.abs(fftpack.ifft(Ar*B))) -> 4
numpy.argmax(numpy.abs(fftpack.ifft(A*Br))) -> 17

再次强调，这两个值对应于你是在解释a的位移还是b的位移。

负共轭是因为卷积会翻转其中一个函数，而在相关性中没有翻转。你可以通过反转其中一个信号然后进行FFT，或者先对信号进行FFT再取负共轭来消除翻转。也就是说，以下是成立的：Ar = -A.conjugate() = fft(a[::-1])