问如何在Numpy中限制互相关窗宽？
EN

Stack Overflow用户

提问于 2015-06-05 23:25:27

回答 8查看 6.8K关注 0票数 7

我正在学习numpy/scipy，来自MATLAB背景。xcorr function in Matlab有一个可选参数" maxlag“，它将lag的范围限制在-maxlag到maxlag之间。如果您正在查看两个非常长的时间序列之间的相互关系，但只对特定时间范围内的相关性感兴趣，这是非常有用的。考虑到交叉相关的计算成本非常高，性能的提升是巨大的。

在numpy/scipy中，似乎有几个计算互相关的选项。numpy.correlate，numpy.convolve，scipy.signal.fftconvolve。如果有人想解释它们之间的区别，我很乐意听到，但主要困扰我的是，它们都没有maxlag特性。这意味着，即使我只想查看滞后在-20000和+20000毫秒之间的两个时间序列之间的相关性，它仍然会计算-100ms和+100毫秒(这是时间序列的长度)之间的每个滞后的相关性。这会带来200倍的性能损失！我是否必须手动重新编码互相关函数才能包含此功能？

numpy

scipy

cross-correlation

python

回答 8

Stack Overflow用户

发布于 2016-01-01 18:47:53

这里有几个函数可以在有限滞后的情况下计算自相关和互相关。选择乘法(在复杂情况下是共轭)的顺序，以匹配numpy.correlate的相应行为。

import numpy as np
from numpy.lib.stride_tricks import as_strided


def _check_arg(x, xname):
    x = np.asarray(x)
    if x.ndim != 1:
        raise ValueError('%s must be one-dimensional.' % xname)
    return x

def autocorrelation(x, maxlag):
    """
    Autocorrelation with a maximum number of lags.

    `x` must be a one-dimensional numpy array.

    This computes the same result as
        numpy.correlate(x, x, mode='full')[len(x)-1:len(x)+maxlag]

    The return value has length maxlag + 1.
    """
    x = _check_arg(x, 'x')
    p = np.pad(x.conj(), maxlag, mode='constant')
    T = as_strided(p[maxlag:], shape=(maxlag+1, len(x) + maxlag),
                   strides=(-p.strides[0], p.strides[0]))
    return T.dot(p[maxlag:].conj())


def crosscorrelation(x, y, maxlag):
    """
    Cross correlation with a maximum number of lags.

    `x` and `y` must be one-dimensional numpy arrays with the same length.

    This computes the same result as
        numpy.correlate(x, y, mode='full')[len(a)-maxlag-1:len(a)+maxlag]

    The return vaue has length 2*maxlag + 1.
    """
    x = _check_arg(x, 'x')
    y = _check_arg(y, 'y')
    py = np.pad(y.conj(), 2*maxlag, mode='constant')
    T = as_strided(py[2*maxlag:], shape=(2*maxlag+1, len(y) + 2*maxlag),
                   strides=(-py.strides[0], py.strides[0]))
    px = np.pad(x, maxlag, mode='constant')
    return T.dot(px)

例如,

In [367]: x = np.array([2, 1.5, 0, 0, -1, 3, 2, -0.5])

In [368]: autocorrelation(x, 3)
Out[368]: array([ 20.5,   5. ,  -3.5,  -1. ])

In [369]: np.correlate(x, x, mode='full')[7:11]
Out[369]: array([ 20.5,   5. ,  -3.5,  -1. ])

In [370]: y = np.arange(8)

In [371]: crosscorrelation(x, y, 3)
Out[371]: array([  5. ,  23.5,  32. ,  21. ,  16. ,  12.5,   9. ])

In [372]: np.correlate(x, y, mode='full')[4:11]
Out[372]: array([  5. ,  23.5,  32. ,  21. ,  16. ,  12.5,   9. ])

(如果numpy本身有这样的特性就好了。)

票数 8

Stack Overflow用户

发布于 2017-12-19 19:20:16

在numpy实现maxlag参数之前，您可以从pycorrelate package使用函数ucorrelate。ucorrelate对numpy数组进行操作，并具有maxlag关键字。它使用for循环实现关联，并使用numba优化执行速度。

具有3个时间延迟的示例-自相关：

import numpy as np
import pycorrelate as pyc

x = np.array([2, 1.5, 0, 0, -1, 3, 2, -0.5])
c = pyc.ucorrelate(x, x, maxlag=3)
c

结果：

Out[1]: array([20,  5, -3])

pycorrelate文档包含显示pycorrelate.ucorrelate和numpy.correlate之间完美匹配的a notebook

票数 5

Stack Overflow用户

发布于 2015-06-06 07:37:23

matplotlib.pyplot提供了类似matlab的语法，用于计算和绘制互相关、自相关等。

您可以使用允许定义maxlags参数的xcorr。

    import matplotlib.pyplot as plt


    import numpy  as np


    data = np.arange(0,2*np.pi,0.01)


    y1 = np.sin(data)


    y2 = np.cos(data)


    coeff = plt.xcorr(y1,y2,maxlags=10)

    print(*coeff)


[-10  -9  -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7
   8   9  10] [ -9.81991753e-02  -8.85505028e-02  -7.88613080e-02  -6.91325329e-02
  -5.93651264e-02  -4.95600447e-02  -3.97182508e-02  -2.98407146e-02
  -1.99284126e-02  -9.98232812e-03  -3.45104289e-06   9.98555430e-03
   1.99417667e-02   2.98641953e-02   3.97518558e-02   4.96037706e-02
   5.94189688e-02   6.91964864e-02   7.89353663e-02   8.86346584e-02
   9.82934198e-02] <matplotlib.collections.LineCollection object at 0x00000000074A9E80> Line2D(_line0)