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

问题描述 投票:0回答:7

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

如何查看是否有时移?如果是的话,多少钱。

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

python numpy signal-processing correlation
7个回答
47
投票

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


15
投票

如果一个被另一个时移,您将看到相关性的峰值。由于计算相关性的成本较高,因此最好使用 FFT。所以,像这样的东西应该有效:

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

time_shift = argmax(abs(c))

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

这不是一种可靠的信号匹配方法,但它快速且简单。


3
投票

这取决于您拥有的信号类型(周期性?...)、两个信号是否具有相同的幅度以及您想要的精度。

highBandWidth 提到的相关函数可能确实适合你。 非常简单,您应该尝试一下。

另一个更精确的选项是我用于高精度谱线拟合的选项:您使用样条曲线对“主”信号进行建模,并用它拟合时移信号(如果需要,同时可能缩放信号)。 这会产生非常精确的时间平移。 这种方法的优点之一是您不必研究相关函数。 例如,您可以使用

interpolate.UnivariateSpline()
(来自 SciPy)轻松创建样条线。 SciPy 返回一个函数,然后可以轻松地使用
optimize.leastsq
() 来拟合该函数。


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

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

请参阅此处的示例


0
投票

我有完全相同的问题,但不喜欢最上面的答案,因为它不会返回关于相关性是否强的信心:如果你最终试图在另一个不包含的数组中找到一个数组它,即使相关性很小,它也会给你答案。

所以,我想出了一个函数,既返回延迟置信度,并将其制作成Python模块

© www.soinside.com 2019 - 2024. All rights reserved.