[python] numpy.correlate를 사용하여 자기 상관을 어떻게 할 수 있습니까?

이 주제에 혼란을 더하는 두 가지가 있다고 생각합니다.

1. 통계 대 신호 처리 정의 : 다른 사람들이 지적했듯이 통계에서 자동 상관 관계를 [-1,1]로 정규화합니다.
2. 부분적 vs 비 부분적 평균 / 분산 : 시계열이 시차> 0에서 이동하면 중복 크기는 항상 <원래 길이보다 작습니다. 원래의 평균과 표준을 사용합니까 (부분이 아님), 아니면 항상 새 평균을 계산하고 항상 변화하는 겹침 (부분)을 사용하여 표준이 차이를 만듭니다. (아마도 이것에 대한 공식적인 용어가 있지만 지금은 “부분”을 사용하겠습니다).

나는 부분적 구별과 부분적 구별이없는 1d 배열의 자동 상관을 계산하는 5 개의 함수를 만들었습니다. 일부는 통계의 공식을 사용하고 일부는 신호 처리 의미에서 상관 관계를 사용하며 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.correlate, FFT)의 계산이 각 중첩에 대해 다른 평균 / 표준을 계산하지 않기 때문입니다.

또한 fft, no padding, non-partial(빨간색 선) 결과는 FFT를 수행하기 전에 시계열을 0으로 채우지 않았으므로 원형 FFT입니다. 왜 다른 곳에서 배웠는지 자세히 설명 할 수 없습니다.

답변

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()

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)

def autocross(x, y):
    c = np.correlate(x, y, "same")
    v = [c[i]/( len(x)-abs( i - (len(x)/2)  ) ) for i in range(len(c))]
    return v