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

Question

我必须比较两个时间 - 电压波形.由于这些波形源的特殊性,其中一个可以是另一个的时移版本.

我怎样才能找到时间转移？如果是的话,它有多少.

我在Python中这样做,并希望使用numpy/scipy库.

Answer 1

scipy提供了一个相关函数,它可以很好地用于小输入,如果你想要非循环相关意味着信号不会环绕.请注意,mode='full'signal.correlation返回的数组的大小是输入信号大小的总和 - 1,因此值来自len(a) + len(b) - 1您所期望的(信号大小-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

两个不同的值对应于换档是否在argmax或a.

如果你想要循环相关和大信号大小,你可以使用卷积/傅立叶变换定理,但需要注意的是相关性与卷积非常相似但不完全相同.

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

这两个值再次对应于您是否解释了转变b或转变a.

负共轭是由于卷积翻转其中一个函数,但在相关性中没有翻转.您可以通过反转其中一个信号然后进行FFT,或者采用信号的FFT然后采用负共轭来撤消翻转.即以下是真的:b

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

对于实值信号,此功能可能更有效.它使用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

它不是一种匹配信号的强大方法,但它快速而简单.