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

find time shift between two similar waveforms

提问by Vishal

I have to compare two time-vs-voltage waveforms. Because of the peculiarity of the sources of these waveforms, one of them can be a time shifted version of the other.

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

How can i find whether there is a time shift? and if yes, how much is it.

我怎样才能知道是否有时移？如果是，是多少。

I am doing this in Python and wish to use numpy/scipy libraries.

我在 Python 中执行此操作并希望使用 numpy/scipy 库。

采纳答案by Gus

scipy provides a correlation function which will work fine for small input and also if you want non-circular correlation meaning that the signal will not wrap around. note that in mode='full', the size of the array returned by signal.correlation is sum of the signal sizes minus one (i.e. len(a) + len(b) - 1), so the value from argmaxis off by (signal size -1 = 20) from what you seem to expect.

scipy 提供了一个相关函数，它适用于小输入，并且如果您想要非圆形相关意味着信号不会环绕。请注意，在 中mode='full'，signal.correlation 返回的数组大小是信号大小的总和减去 1（即len(a) + len(b) - 1），因此from 的值与argmax您似乎期望的值相差(signal size -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

The two different values correspond to whether the shift is in aor b.

两个不同的值对应于移位是 ina还是b。

If you want circular correlation and for big signal size, you can use the convolution/Fourier transform theorem with the caveat that correlation is very similar to but not identical to convolution.

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

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

again the two values correspond to whether your interpreting a shift in aor a shift in b.

再次将两个值符合您是否解释在移位a或移位b。

The negative conjugation is due to convolution flipping one of the functions, but in correlation there is no flipping. You can undo the flipping by either reversing one of the signals and then taking the FFT, or taking the FFT of the signal and then taking the negative conjugate. i.e. the following is true: Ar = -A.conjugate() = fft(a[::-1])

负共轭是由于卷积翻转了函数之一，但相关地没有翻转。您可以通过反转其中一个信号然后采用 FFT，或采用信号的 FFT 然后采用负共轭来撤消翻转。即以下是正确的：Ar = -A.conjugate() = fft(a[::-1])

回答by highBandWidth

If one is time-shifted by the other, you will see a peak in the correlation. Since calculating the correlation is expensive, it is better to use FFT. So, something like this should work:

如果一个被另一个时移，您将看到相关性的峰值。由于计算相关性很昂贵，因此最好使用 FFT。所以，这样的事情应该工作：

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

time_shift = argmax(abs(c))

回答by Eric O Lebigot

It depends on the kind of signal you have (periodic?…), on whether both signals have the same amplitude, and on what precision you are looking for.

这取决于您拥有的信号类型（周期性？...），取决于两个信号是否具有相同的幅度，以及您正在寻找的精度。

The correlation function mentioned by highBandWidth might indeed work for you. It is simple enough that you should give it a try.

highBandWidth 提到的相关函数可能确实适合您。这很简单，你应该试一试。

Another, more precise option is the one I use for high-precision spectral line fitting: you model your "master" signal with a spline and fit the time-shifted signal with it (while possibly scaling the signal, if need be). This yields very precise time shifts. One advantage of this approach is that you do not have to study the correlation function. You can for instance create the spline easily with interpolate.UnivariateSpline()(from SciPy). SciPy returns a function, which is then easily fitted with optimize.leastsq().

另一个更精确的选项是我用于高精度谱线拟合的选项：您使用样条对“主”信号进行建模，并用它拟合时移信号（同时可能缩放信号，如果需要）。这会产生非常精确的时间偏移。这种方法的一个优点是您不必研究相关函数。例如，您可以使用interpolate.UnivariateSpline()（来自 SciPy）轻松创建样条曲线。SciPy 返回一个函数，然后很容易用optimize.leastsq()拟合它。

回答by Eryk Sun

This function is probably more efficient for real-valued signals. It uses rfft and zero pads the inputs to a power of 2 large enough to ensure linear (i.e. non-circular) correlation:

对于实值信号，此函数可能更有效。它使用 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

The return value is length M = len(x) + len(y) - 1(hacked together with hstackto remove the extra zeros from rounding up to a power of 2). The non-negative lags are cxy[0], cxy[1], ..., cxy[len(x)-1], while the negative lags are cxy[-1], cxy[-2], ..., cxy[-len(y)+1].

返回值是长度M = len(x) + len(y) - 1（与 hack 一起hstack删除额外的零，从四舍五入到 2 的幂）。非负滞后为cxy[0], cxy[1], ..., cxy[len(x)-1]，而负滞后为cxy[-1], cxy[-2], ..., cxy[-len(y)+1]。

To match a reference signal, I'd compute rfft_xcorr(x, ref)and look for the peak. For example:

为了匹配参考信号，我会计算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

It's not a robust way to match signals, but it is quick and easy.

这不是匹配信号的可靠方法，但它快速而简单。

回答by FFT

Here's another option:

这是另一种选择：

from scipy import signal, fftpack

def get_max_correlation(original, match):
    z = signal.fftconvolve(original, match[::-1])
    lags = np.arange(z.size) - (match.size - 1)
    return ( lags[np.argmax(np.abs(z))] )

回答by Mauricio Redaelli

Blockquote

块引用

(A very late answer) to find the time-shift between two signals: use the time-shift property of FTs, so the shifts can be shorter than the sample spacing, then compute the quadratic difference between a time-shifted waveform and the reference waveform. It can be useful when you have n shifted waveforms with a multiplicity in the shifts, like n receivers equally spaced for a same incoming wave. You can also correct dispersion substituting a static time-shift by a function of frequency.

（一个很晚的答案）找到两个信号之间的时移：使用 FT 的时移特性，因此移位可以比样本间隔短，然后计算时移波形与参考之间的二次差波形。当您有 n 个具有多个移位的移位波形时，它会很有用，例如 n 个接收器等距接收相同的传入波。您还可以通过频率函数代替静态时移来校正色散。

The code goes like this:

代码是这样的：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftshift, fftfreq
from scipy import signal

#  generating a test signal
dt = 0.01
t0 = 0.025
n = 512
freq = fftfreq(n, dt)

time = np.linspace(-n * dt / 2, n * dt / 2, n)
y = signal.gausspulse(time, fc=10, bw=0.3) + np.random.normal(0, 1, n) / 100
Y = fft(y)
# time-shift of 0.235; could be a dispersion curve, so y2 would be dispersive
Y2 = Y * np.exp(-1j * 2 * np.pi * freq * 0.235)  
y2 = ifft(Y2).real

# scan possible time-shifts
error = []
timeshifts = np.arange(-100, 100) * dt / 2  # could be dispersion curves instead
for ts in timeshifts:
    Y2_shifted = Y2 * np.exp(1j * 2 * np.pi * freq * ts)
    y2_shifted = ifft(Y2_shifted).real
    error.append(np.sum((y2_shifted - y) ** 2))

# show the results
ts_final = timeshifts[np.argmin(error)]
print(ts_final)

Y2_shifted = Y2 * np.exp(1j * 2 * np.pi * freq * ts_final)
y2_shifted = ifft(Y2_shifted).real

plt.subplot(221)
plt.plot(time, y, label="y")
plt.plot(time, y2, label="y2")
plt.xlabel("time")
plt.legend()

plt.subplot(223)
plt.plot(time, y, label="y")
plt.plot(time, y2_shifted, label="y_shifted")
plt.xlabel("time")
plt.legend()

plt.subplot(122)
plt.plot(timeshifts, error, label="error")
plt.xlabel("timeshifts")
plt.legend()

plt.show()

See an example here

在此处查看示例