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

Posted 2023-02-14

【中文标题】找出两个相似波形之间的时间偏移

【英文标题】：find time shift between two similar waveforms

【发布时间】：2011-06-08 23:47:34

【问题描述】：

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

如何确定是否存在时移？如果是的话，多少钱。

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

【参考方案1】：

scipy 提供了一个相关函数，该函数适用于小输入，并且如果您想要非循环相关，这意味着信号不会环绕。请注意，在mode='full' 中，signal.correlation 返回的数组大小是信号大小的总和减一（即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

这两个不同的值对应shift是在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])

【讨论】：

感谢您的回答。这是我第一次看到有意义的东西。现在还有一个问题，取决于时移值的“符号”，我将减去或添加时移。如何获得标志？

等等...你为什么需要负片？我认为你不需要负面的。令 x(t) 具有变换 X(f)。通过时间反转，x(-t) 具有变换 X(-f)。如果 x(t) 是实数，则 X(-f) = conj(X(f))。因此，如果 x(t) 是实数，则 x(-t) 具有变换 conj(X(f))。没有负面的。

@Steve：谢谢。我昨晚推导它时犯了一个错误。

感谢您的回答 - 它也帮助我解决了我的问题。

@SteveTjoa Vishal 注意到的是 signal.correlate 不假设信号是周期性的，因此返回正向或负向偏移，而第二种方法总是返回正向偏移，这是可以的，因为信号是应该是周期性的。

【参考方案2】：

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

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

time_shift = argmax(abs(c))

【讨论】：

我尝试按照您的建议进行操作，因为手头的案例给出了错误的结果。示例： >>> a21 数组（[0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0 ]) >>> a22 数组（[0, 0, 0, 0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0 ]) >>> fa21 = np.fft.fft(a21) >>> fa22 = np.fft.fft(a22) >>> c = np.fft.ifft(fa21 * fa22) >>> time_shift = np. argmax(abs(c)) >>> time_shift 20 正如你所看到的，实际的时间偏移是 4 点而不是 20。我在这里遗漏了什么吗？

-1。不正确，因为 c 只是 a 与 b 卷积，不相关。时间倒转会使事情变得一团糟，不会产生预期的结果。

你说得对，史蒂夫。我把答案写成一个粗略的想法。我已经更正了它以反映共轭。

感谢您的编辑。 （这仅适用于真实信号，但我想我们可以假设。）

有没有办法找出领先的信号？

【参考方案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

这不是一种可靠的信号匹配方式，但它既快速又简单。

【参考方案4】：

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

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

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

【讨论】：

谢谢！我刚刚使用了 optimize.leastsq：我不知道这对于时移来说是易于处理的；比卷积方法容易得多。你知道是否有任何关于 optimize.leastsq 工作原理的参考资料吗？我认为最小二乘必须处理输入基函数的线性组合。

在documentation 中写道“leastsq”是MINPACK 的lmdif 和lmder 算法的封装。”您可以在MINPACK 的代码中找到更多信息：netlib.org/minpack/lmdif.f 和netlib.org/minpack/lmder.f。

【参考方案5】：

这是另一个选择：

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))] )

【讨论】：

工作，但 seems completely equivalent 到 scipy.signal.correlate() 从 Gus answer 默认使用 scipy.signal.fftconvolve 尽快（即一旦二次次伤害很快)。

当数据为例如增加。 a= [ 0 2 4 6 8 8 8 8 8 10 12 14 16 16 16 16 16 17 18 19 20] b=[-4 -3 -2 -1 0 2 4 6 8 8 8 8 8 10 12 14 16 16 16 16 16] get_max_correlation(a,b) -> 0，当应用 a=numpy.gradient(a) b=numpy.gradient(b) 它正确返回 get_max_correlation(a,b) -> -4

【参考方案6】：

块引用

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

代码如下：

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