Cython与numpy速度

Question

我在我的Python程序中使用Cython来计算相关性。我有两个音频数据集，我需要知道它们之间的时间差。第二个数据集是根据开始时间切割的，然后在第一个数据集上滑动。这里有两个循环：一个是滑动数据集，另一个是在那个点计算相关性。这个方法效果很好，准确度也足够。

问题是，使用纯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

Answer 1

这种情况的关键是找到一种“分而治之”的方法。

现在，你是在每个位置滑动并检查每个点，这样做的效率很低，实际上是 O( n ^ 2 ) 的操作。

你需要减少对每个点的检查和每个位置的比较，把工作量降到最低，以便更快地判断是否不匹配。

比如，你可以设置一个较短的“这是否接近？”的过滤器，先检查前几个位置。如果相关性超过某个阈值，就继续检查，否则就放弃，换个方向。

你还可以选择“每隔8个位置检查一次”，然后乘以8。如果这个值太低，就跳过，继续下一个。如果这个值足够高，就检查所有的值，看看是否找到了最大值。

问题在于进行这些乘法运算所需的时间——（f[<unsigned int>(i+j)] * g[j]）。实际上，你是在填充一个大矩阵，计算所有这些乘积，然后选择和最大的那一行。你并不想计算“所有”的乘积，而只需要计算足够的乘积，以确保找到了最大和。

寻找最大值的问题在于，你必须把所有的值相加，才能判断哪个最大。如果你能把这个问题转化为最小化问题，那就更简单了。一旦中间结果超过某个阈值，就可以放弃计算乘积和和。

（我觉得这个方法可能有效，但我还没试过。）

如果你用 max(g)-g[j] 来处理负数，你就要找最小值，而不是最大值。你可以先计算第一个位置的相关性。任何和更大的值都可以立即停止计算——不再进行乘法或加法，直接换个位置。

Answer 2

你可以把外层循环中的范围（size2）提取出来。

你可以用sum()这个函数来计算当前的相关性，而不需要用循环。

你可以把相关性和延迟存储在一个列表里，然后用max()函数来找出最大的那个。

Answer 3

编辑：
现在有了一个新的方法 scipy.signal.fftconvolve，这是我下面描述的基于FFT的卷积方法的推荐做法。我会保留原来的回答来解释速度问题，但在实际使用中，建议使用 scipy.signal.fftconvolve。

原回答：
使用快速傅里叶变换（FFTs）和卷积定理，可以大幅提高速度，把问题的复杂度从 O(n^2) 降到 O(n log n)。这对于像你这样长的数据集特别有用，速度提升可以达到上千倍甚至更多，具体取决于数据的长度。这个方法也很简单：只需对两个信号进行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则把波形视为周期性的。如果这个问题影响到结果，可以通过零填充来解决。