如何限制Numpy中的互相關窗口寬度？

Question

我正在學習numpy/scipy，來自MATLAB背景。 xcorr function in Matlab有一個可選參數「maxlag」，它限制從-maxlag到maxlag的滯後範圍。如果您正在查看兩個很長時間序列之間的互相關，但只對某個時間範圍內的相關性感興趣，這非常有用。考慮到互相關計算的成本非常昂貴，性能的提高是巨大的。

在numpy/scipy中，似乎有幾種計算互相關的選項。 numpy.correlate，numpy.convolve，scipy.signal.fftconvolve。如果有人想解釋它們之間的差異，我很樂意聽到，但主要是令我感到困擾的是，他們都沒有maxlag功能。這意味着，即使我只想看到兩個時間序列之間的相關性，例如，滯後時間在-100到+100 ms之間，它仍然會計算-20000到+20000 ms之間的每個滯後的相關性（這是時間序列）。這給了200倍的性能打擊！我是否必須手動重新編碼互相關函數才能包含此功能？

Answer 1

matplotlib.pyplot提供MATLAB像語法computating和互相關的繪圖，自相關等

可以使用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) 

Answer 2

我想我已經找到了解決辦法，因爲我面臨同樣的問題：

如果你有兩個向量x和任何長度爲N y，並且想使用的窗口互相關固定len個m，你可以這樣做：

x = <some_data> 
y = <some_data> 

# Trim your variables 
x_short = x[window:] 
y_short = y[window:] 

# do two xcorrelations, lagging x and y respectively 
left_xcorr = np.correlate(x, y_short) #defaults to 'valid' 
right_xcorr = np.correlate(x_short, y) #defaults to 'valid' 

# combine the xcorrelations 
# note the first value of right_xcorr is the same as the last of left_xcorr 
xcorr = np.concatenate(left_xcorr, right_xcorr[1:]) 

記住，如果你想有一個有限相關，你可能需要normalise變量

Answer 3

這裏是另一個答案，從here來源，似乎在緣速度比np.correlate並有回國標準化相關的好處：

def rolling_window(self, a, window): 
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) 
    strides = a.strides + (a.strides[-1],) 
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) 

def xcorr(self, x,y): 

    N=len(x) 
    M=len(y) 
    meany=np.mean(y) 
    stdy=np.std(np.asarray(y)) 
    tmp=self.rolling_window(np.asarray(x),M) 
    c=np.sum((y-meany)*(tmp-np.reshape(np.mean(tmp,-1),(N-M+1,1))),-1)/(M*np.std(tmp,-1)*stdy) 

    return c   

Answer 4

這裏有幾個功能，自動計算，並與有限的互相關滯後。選擇乘法（和複合情況下的共軛）的順序以匹配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 5

直到numpy的實現maxlag參數，您可以從使用功能ucorrelatepycorrelate package。 ucorrelate在numpy數組上運行，並有一個maxlag關鍵字。它通過使用for-loop實現相關性，並使用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文檔包含a notebook表示pycorrelate.ucorrelate和numpy.correlate之間完美匹配：

Answer 6

正如我在這裏回答的那樣，https://stackoverflow.com/a/47897581/5122657 matplotlib.xcorr有maxlags參數。它實際上是numpy.correlate的包裝，因此不會節省性能。不過，它給出了Matlab的互相關函數給出的完全相同的結果。下面我編輯了matplotlib中的代碼，以便它只返回相關性。原因是，如果我們按原樣使用matplotlib.corr，它也會返回圖表。問題是，如果我們把複雜的數據類型作爲參數放入它，當matplotlib試圖繪製圖時，我們會得到「複雜的實體數據類型」警告。

<!-- language: python --> 

import numpy as np 
import matplotlib.pyplot as plt 

def xcorr(x, y, maxlags=10): 
    Nx = len(x) 
    if Nx != len(y): 
     raise ValueError('x and y must be equal length') 

    c = np.correlate(x, y, mode=2) 

    if maxlags is None: 
     maxlags = Nx - 1 

    if maxlags >= Nx or maxlags < 1: 
     raise ValueError('maxlags must be None or strictly positive < %d' % Nx) 

    c = c[Nx - 1 - maxlags:Nx + maxlags] 

    return c