16

我在我的 python 程序中使用 cython 进行相关性计算。我有两个音频数据集，我需要知道它们之间的时间差。第二组根据开始时间进行切割，然后滑过第一组。有两个 for 循环：一个滑动集合，内部循环计算该点的相关性。这种方法效果很好，而且足够准确。

问题是，对于纯 python，这需要超过一分钟。使用我的 cython 代码，大约需要 17 秒。这还是太多了。您是否有任何提示如何加速此代码：

import numpy as np
cimport numpy as np

cimport cython

FTYPE = np.float
ctypedef np.float_t FTYPE_t

@cython.boundscheck(False)
def delay(np.ndarray[FTYPE_t, ndim=1] f, np.ndarray[FTYPE_t, ndim=1] g):
    cdef int size1 = f.shape[0]
    cdef int size2 = g.shape[0]
    cdef int max_correlation = 0
    cdef int delay = 0
    cdef int current_correlation, i, j

    # Move second data set frame by frame
    for i in range(0, size1 - size2):
        current_correlation = 0

        # Calculate correlation at that point
        for j in range(size2):
            current_correlation += f[<unsigned int>(i+j)] * g[j]

        # Check if current correlation is highest so far
        if current_correlation > max_correlation:
            max_correlation = current_correlation
            delay = i

    return delay
4

3 回答 3

37

编辑：
现在scipy.signal.fftconvolve这将是我在下面描述的基于 FFT 的卷积方法的首选方法。我将留下原始答案来解释速度问题，但在实践中使用scipy.signal.fftconvolve.

原始答案：
使用FFT和卷积定理将问题从 O(n^2) 转换为 O(n log n)，从而显着提高速度。这对于像您这样的长数据集特别有用，并且可以提供 1000 秒或更多的速度增益，具体取决于长度。这也很容易做到：只需对两个信号进行 FFT、乘法和逆 FFT 乘积。numpy.correlate在互相关例程中不使用 FFT 方法，最好与非常小的内核一起使用。

这是一个例子

from timeit import Timer
from numpy import *

times = arange(0, 100, .001)

xdata = 1.*sin(2*pi*1.*times) + .5*sin(2*pi*1.1*times + 1.)
ydata = .5*sin(2*pi*1.1*times)

def xcorr(x, y):
    return correlate(x, y, mode='same')

def fftxcorr(x, y):
    fx, fy = fft.fft(x), fft.fft(y[::-1])
    fxfy = fx*fy
    xy = fft.ifft(fxfy)
    return xy

if __name__ == "__main__":
    N = 10
    t = Timer("xcorr(xdata, ydata)", "from __main__ import xcorr, xdata, ydata")
    print 'xcorr', t.timeit(number=N)/N
    t = Timer("fftxcorr(xdata, ydata)", "from __main__ import fftxcorr, xdata, ydata")
    print 'fftxcorr', t.timeit(number=N)/N

它给出了每个周期的运行时间（以秒为单位，对于 10,000 长波形）

xcorr 34.3761689901
fftxcorr 0.0768054962158

很明显 fftxcorr 方法要快得多。

如果您绘制结果，您会发现它们在接近零时移时非常相似。但是请注意，随着您离得越远，xcorr 会减少，而 fftxcorr 不会。这是因为如何处理波形移动时不重叠的波形部分有点模棱两可。xcorr 将其视为零，FFT 将波形视为周期性的，但如果这是一个问题，则可以通过零填充来修复。

于 2009-08-04T17:40:00.020 回答
2

这种事情的诀窍是找到一种分而治之的方法。

目前，您正在滑动到每个位置并检查每个位置的每个点——实际上是O ( n ^ 2 ) 操作。

您需要减少对每个点的检查和每个位置的比较，以减少确定不匹配的工作。

例如，您可以使用更短的“这是否更接近？” 检查前几个位置的过滤器。如果相关性高于某个阈值，则继续前进，否则放弃并继续前进。

你可以有一个“每第 8 个位置检查一次”，然后乘以 8。如果这太低，请跳过它并继续。如果这足够高，则检查所有值以查看是否找到最大值。

问题是进行所有这些乘法所需的时间—— ( f[<unsigned int>(i+j)] * g[j]) 实际上，您正在用所有这些乘积填充一个大矩阵并选择总和最大的行。您不想计算“所有”产品。足够的产品可以确保您找到最大金额。

找到最大值的问题是你必须把所有的东西都加起来看看它是否最大。如果您可以将其转化为最小化问题，那么一旦中间结果超过阈值，就更容易放弃计算产品和求和。

（我认为这可能有效。我还没有尝试过。）

如果您曾经使用max(g)-g[j]负数，您会寻找最小的，而不是最大的。您可以计算第一个位置的相关性。任何总和更大的值都可以立即停止——不再为那个偏移量乘法或加法，转移到另一个。

于 2009-07-29T12:58:22.203 回答
2
  • 您可以从外部循环中提取 range(size2)
  • 您可以使用 sum() 而不是循环来计算 current_correlation
  • 您可以将相关性和延迟存储在列表中，然后使用 max() 获取最大的
于 2009-07-29T13:04:13.700 回答