如何使用numpy.correlate进行自相关？

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]

当然,您需要进行错误检查以确保它x实际上是一维数组.此外,这种解释可能不是数学上最严格的.由于卷积的定义使用了无限性,我一直在无穷无尽,但这并不一定适用于自相关.因此,这种解释的理论部分可能略显不足,但希望实际结果是有帮助的.这些 关于自相关的页面非常有用,如果您不介意使用符号和繁重的概念,可以为您提供更好的理论背景.

自相关有两个版本:统计和卷积.除了一些细节之外,它们都是相同的:统计版本被标准化为区间[-1,1].以下是如何进行统计分析的示例:

def acf(x, length=20):
    return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:])[0,1]  \
        for i in range(1, length)])

使用numpy.corrcoef函数代替numpy.correlate计算滞后t的统计相关性:

def autocorr(x, t=1):
    return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))

由于我遇到了同样的问题,我想与您分享几行代码.事实上,目前有一些关于stackoverflow中自相关的相似帖子.如果将自相关定义为a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2)[这是IDL的a_correlate函数中给出的定义,并且它与我在问题#12269834的答案2中看到的内容一致],则以下似乎给出了正确的结果:

import numpy as np
import matplotlib.pyplot as plt

# generate some data
x = np.arange(0.,6.12,0.01)
y = np.sin(x)
# y = np.random.uniform(size=300)
yunbiased = y-np.mean(y)
ynorm = np.sum(yunbiased**2)
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm
# use only second half
acor = acor[len(acor)/2:]

plt.plot(acor)
plt.show()

如你所见,我用sin曲线和均匀的随机分布测试了这个,两个结果看起来都像我期望的那样.请注意,我使用mode="same"而不是mode="full"像其他人那样.

我认为有两件事会增加这个话题的混乱:

统计与信号处理定义:正如其他人所指出的,在统计中我们将自相关归一化为[-1,1].

部分与非部分均值/方差:当时间序列在滞后> 0时移位时,它们的重叠大小将始终<原始长度.我们是使用原始(非部分)的均值和标准,还是总是使用不断变化的重叠(部分)来计算新均值和标准值.(这可能是一个正式的术语,但我现在要使用"部分").

我创建了5个函数来计算1d数组的自相关,部分与非部分区别.有些使用统计公式,有些使用信号处理意义上的相关,也可以通过FFT完成.但是所有结果都是统计定义中的自相关,因此它们说明了它们如何相互关联.代码如下:

import numpy
import matplotlib.pyplot as plt

def autocorr1(x,lags):
    '''numpy.corrcoef, partial'''

    corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags]
    return numpy.array(corr)

def autocorr2(x,lags):
    '''manualy compute, non partial'''

    mean=numpy.mean(x)
    var=numpy.var(x)
    xp=x-mean
    corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags]

    return numpy.array(corr)

def autocorr3(x,lags):
    '''fft, pad 0s, non partial'''

    n=len(x)
    # pad 0s to 2n-1
    ext_size=2*n-1
    # nearest power of 2
    fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int')

    xp=x-numpy.mean(x)
    var=numpy.var(x)

    # do fft and ifft
    cf=numpy.fft.fft(xp,fsize)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real
    corr=corr/var/n

    return corr[:len(lags)]

def autocorr4(x,lags):
    '''fft, don't pad 0s, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean

    cf=numpy.fft.fft(xp)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real/var/len(x)

    return corr[:len(lags)]

def autocorr5(x,lags):
    '''numpy.correlate, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean
    corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)

    return corr[:len(lags)]


if __name__=='__main__':

    y=[28,28,26,19,16,24,26,24,24,29,29,27,31,26,38,23,13,14,28,19,19,\
            17,22,2,4,5,7,8,14,14,23]
    y=numpy.array(y).astype('float')

    lags=range(15)
    fig,ax=plt.subplots()

    for funcii, labelii in zip([autocorr1, autocorr2, autocorr3, autocorr4,
        autocorr5], ['np.corrcoef, partial', 'manual, non-partial',
            'fft, pad 0s, non-partial', 'fft, no padding, non-partial',
            'np.correlate, non-partial']):

        cii=funcii(y,lags)
        print(labelii)
        print(cii)
        ax.plot(lags,cii,label=labelii)

    ax.set_xlabel('lag')
    ax.set_ylabel('correlation coefficient')
    ax.legend()
    plt.show()

这是输出数字:

我们没有看到所有5条线,因为它们中的3条重叠(在紫色处).重叠都是非部分自相关.这是因为来自信号处理方法(np.correlateFFT)的计算不会为每个重叠计算不同的均值/标准.

另请注意,fft, no padding, non-partial(红线)结果不同,因为它在进行FFT之前没有用0填充时间序列,因此它是循环FFT.我无法详细解释为什么,这就是我从其他地方学到的东西.

您的问题1已在此处的几个优秀答案中进行了广泛讨论.

我想与您分享几行代码,这些代码允许您仅根据自相关的数学属性计算信号的自相关性.也就是说,可以通过以下方式计算自相关:

从信号中减去均值并获得无偏信号

计算无偏信号的傅里叶变换

通过取无偏信号的傅立叶变换的每个值的平方范数来计算信号的功率谱密度

计算功率谱密度的逆傅里叶变换

将功率谱密度的傅立叶逆变换归一化无偏信号的平方和,并仅取得结果矢量的一半

执行此操作的代码如下:

def autocorrelation (x) :
    """
    Compute the autocorrelation of the signal, based on the properties of the
    power spectral density of the signal.
    """
    xp = x-np.mean(x)
    f = np.fft.fft(xp)
    p = np.array([np.real(v)**2+np.imag(v)**2 for v in f])
    pi = np.fft.ifft(p)
    return np.real(pi)[:x.size/2]/np.sum(xp**2)