避免for循环的多个序列的互相关

Question

阅读 this 并尝试了 np.correlate 和 cv2.matchTemplate 我仍然有一个我似乎无法解决的问题。

我有两个 numpy 数组，每个数组的形状都为 (6000,50)。 6000 个序列，每 50 个值。现在我想做这个数组的两个一维序列的互相关来检测时移。我简单地尝试了openCV，但对我来说这会返回一个数字（我希望相关性最高），所以现在我像这样使用numpy.correlate：

np.correlate(x[2500], y[2500], mode='same')

（在互相关图中，我不是在寻找最高峰，而是在寻找使用this 的第一个峰。请参见图中的示例）

如您所料，我想对所有 6000 个序列执行此操作，但希望避免迭代。我希望这会奏效：

np.correlate(x, y, mode='same')

但这给了我以下错误：ValueError: object too deep for desired array。

NumPy 或 OpenCV 是否有任何改变。还是我必须这样做:(

for i in range(x.shape[0]):
    np.correlate(x[i], y[i], mode='same')

Answer 1

scipy.ndimage.correlate1d 似乎是你所追求的，但它只在第一个数组上广播，第二个必须是严格的一维，所以那里没有运气。而scipy.signal 中的函数是多维相关的，而不是像你所追求的一维。因此，堆栈中似乎没有任何东西可以解决您的问题。

只是为了好玩，您总是可以使用 FFT 和 cross-correlation theorem：

def correlate1(a, b):
    c = np.empty_like(a)
    for j in range(len(a)):
        c[j] = np.correlate(a[j], b[j], 'same')
    return c

def correlate2(a, b):
    n = a.shape[-1]
    a_fft = np.fft.fft(a, n=2*n)
    b_fft = np.fft.fft(b, n=2*n)
    cc = np.fft.ifft(a_fft * b_fft.conj()).real
    return np.concatenate((cc[..., -n//2:], cc[..., :(n-1)//2 + 1]), axis=-1)

对于您的用例，这不是一个好主意：

In [11]: a = np.random.rand(6000, 50)
    ...: b = np.random.rand(6000, 50)
    ...: 

In [12]: np.allclose(correlate1(a, b), correlate2(a, b))
Out[12]: True

In [13]: %timeit correlate1(a, b)
10 loops, best of 3: 37.5 ms per loop

In [14]: %timeit correlate2(a, b)
10 loops, best of 3: 71.8 ms per loop

但该方法确实有其优点，主要用于较大的序列：

In [15]: a = np.random.rand(50, 6000)
    ...: b = np.random.rand(50, 6000)
    ...: 

In [16]: %timeit correlate1(a, b)
1 loops, best of 3: 516 ms per loop

In [17]: %timeit correlate2(a, b)
10 loops, best of 3: 89.2 ms per loop