找出两个类似波形之间的时间偏移量。

Question

找出两个类似波形之间的时间偏移量。

pythonnumpysignal-processingcorrelation

28

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

如何找出是否存在时间偏移？如果有，它是多少。

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

- Vishal

6个回答

15

如果其中一个序列被另一个序列进行了时间平移，你会看到它们之间的相关性会有一个峰值。由于计算相关性是很昂贵的，所以最好使用快速傅里叶变换（FFT）来计算。因此，以下代码可能会起作用：

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

time_shift = argmax(abs(c))

- highBandWidth

1

我尝试了你建议的方法，但对于这种情况它给出了错误的结果。例如：

>>> a21 
array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0,   0])
>>> a22
array([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个点。我有什么遗漏吗？ - Vishal

1

-1. 这是错误的，因为 c 只是 a 和 b 卷积，而不是相关。时间反转会弄乱事情，而且不会得到想要的结果。 - Steve Tjoa

1

你说得对，Steve。我之前写的答案只是一个大概的想法。我已经更正了它以反映正确的动词变化形式。 - highBandWidth

2

有没有办法找出哪个信号是领先的？ - Shashank Sawant

我认为如果实际Ifft输出的峰值在开头，信号1相对于信号2会发生移位。如果峰值在输出末尾，则反之。 - thomas.cloud

显示剩余2条评论

9

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

这不是一种强大的信号匹配方式，但它快速简便。

- Eryk Sun

3

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

highBandWidth提到的相关函数可能适合您。它足够简单，您应该尝试一下。

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

- Eric O Lebigot

谢谢！我刚刚使用了optimize.leastsq：我不知道这对于时间偏移是可行的；比卷积方法容易得多。您是否知道有关optimize.leastsq如何工作的任何参考资料？我认为最小二乘法必须使用输入基函数的线性组合。 - Jason S

1

在文档中，我们可以看到“leastsq”是MINPACK的lmdif和lmder算法的包装器。您可以在MINPACK的代码中找到更多信息：http://www.netlib.org/minpack/lmdif.f和http://www.netlib.org/minpack/lmder.f。 - Eric O. Lebigot

这真的很有趣，但我正在努力编写它。我找到的所有optimize.leastsq()示例似乎都是关于从适合于函数的实际数据提取模型系数。我已经从“主信号”中得到了样条函数，但如何将另一个数据集拟合到该函数并提取时间偏移量呢？ - vantom

3

这里还有另一个选择：

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

- FFT

能够正常工作，但是似乎与Gus的回答中提到的scipy.signal.correlate()函数完全等价。该函数默认使用scipy.signal.fftconvolve，因为它更快（即在二次时间受损之前就更快）。 - Gabriel Devillers

当数据例如递增时，它会与Gus的答案一样失败。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。 - Marco

1

引用

（一个非常晚的答案）要找到两个信号之间的时间偏移量：使用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()

在此处查看示例

- Mauricio Redaelli

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Gus · Accepted Answer

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

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