如何限制Numpy中的互相关窗口宽度？

Question

我正在学习numpy/scipy,来自MATLAB背景.Matlab中的xcorr函数有一个可选参数"maxlag",它将滞后范围限制在-maxlag到maxlag之间.如果您正在查看两个非常长的时间序列之间的互相关,但仅对某个时间范围内的相关性感兴趣,这非常有用.考虑到互相关的计算成本非常高,性能的提升是巨大的.

在numpy/scipy中,似乎有几种计算互相关的选项. numpy.correlate,numpy.convolve,scipy.signal.fftconvolve.如果有人想解释它们之间的区别,我会很高兴听到,但主要是令我不安的是它们都没有maxlag功能.这意味着即使我只想看到两个时间序列之间的相关性,例如在-100和+100毫秒之间,它仍将计算-20000和+20000毫秒之间的每个滞后的相关性(这是时间序列).这样可以获得200倍的性能!我是否必须手动重新编码互相关函数以包含此功能？

Answer 1

以下是一些用于计算有限滞后的自相关和互相关的函数.选择乘法(和复杂情况下的共轭)的顺序以匹配相应的行为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本身拥有这样的功能会很高兴.)

Answer 2

直到numpy的实现maxlag参数，你可以使用函数ucorrelate从pycorrelate包。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：

Answer 3

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

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

    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)