我如何使用numpy.correlate做自相关?
我需要做一组数字的自动关联,据我所知,它只是集合与自身的相关性。
我试过用numpy的相关函数,但我不相信结果,因为它几乎总是给出一个向量,其中第一个数不是最大的,它应该是这样的。
所以,这个问题真的是两个问题:
- numpy.correlate究竟是干什么的?
- 我怎样才能使用它(或其他)做自相关?
为了回答你的第一个问题, numpy.correlate(a, v, mode)
正在执行a
与v
的相反的卷积,并给出由指定模式截取的结果。 卷积的定义 C(t)=Σ- ∞<i <∞a i v t + i其中-∞<t <∞允许从-∞到∞的结果,但显然不能存储无限长arrays。 所以它必须被剪裁,这就是模式进来的地方。有3种不同的模式:完整,相同,有效:
- “全”模式返回结果为每个
t
,其中a
和v
有一些重叠。 - “相同”模式返回与最短向量(
a
或v
)相同长度的结果。 - 只有当
a
和v
完全重叠时,“有效”模式才返回结果。numpy.convolve
的文档给出了关于模式的更多细节。
对于你的第二个问题,我认为numpy.correlate
给你的自相关,它只是给你一点点。 自相关用于发现在某个时间差异下信号或函数与自身有多相似。 在时间差为0时,自相关应该是最高的,因为信号与自身相同,所以您预期自相关结果数组中的第一个元素将是最大的。 然而,相关性不是从0的时间差开始的,它始于负的时间差,接近0,然后变为正值。 那就是,你在期待:
自相关(a)=Σ- ∞<i <∞a i v t + i其中0 <= t <∞
但是你得到的是:
自相关(a)=Σ- ∞<i <∞a i v t + i其中-∞<t <∞
你需要做的是把你的相关结果的后半部分,这应该是你正在寻找的自相关。 一个简单的python函数可以做到这一点:
def autocorr(x): result = numpy.correlate(x, x, mode='full') return result[result.size/2:]
当然,您将需要错误检查来确保x
实际上是一个一维数组。 而且,这个解释可能不是最严谨的math。 因为卷积的定义使用了它,所以我一直在讨论无穷大,但这并不一定适用于自相关。 所以,这个解释的理论部分可能会有些诡异,但希望实际结果是有帮助的。 这些关于自相关的页面是非常有帮助的,如果你不介意使用符号和沉重的概念,可以给你一个更好的理论背景。
使用numpy.corrcoef
函数代替numpy.correlate
来计算滞后时间t的统计相关性:
def autocorr(x, t=1): numpy.corrcoef(numpy.array([x[0:len(x)-t], x[t:len(x)]]))
自相关有两个版本:统计和卷积。 除了一些细节之外,它们都是一样的:前者被归一化为[-1,1]区间。 这是一个如何做统计的例子:
def acf(x, length=20): return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:]) \ for i in range(1, length)])
当我遇到同样的问题时,我想与您分享几行代码。 事实上,现在有几个相当类似的post关于自动相关在stackoverflow中。 如果将自相关定义为a(x, L) = sum(k=0,NL-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()
正如你所看到的,我已经用一条正弦曲线和一个统一的随机分布来testing这个结果,而且这两个结果看起来像我期望的那样。 请注意,我使用mode="same"
而不是mode="full"
。
您的问题1已经在这里的几个优秀的答案已经广泛讨论。
我想与你分享几行代码,这些代码可以让你只根据自相关的math特性来计算信号的自相关。 也就是说,自相关可以通过以下方式计算:
-
从信号中减去平均值并获得无偏信号
-
计算无偏信号的傅里叶变换
-
通过取无偏信号的傅立叶变换的每个值的平方范数来计算信号的功率谱密度
-
计算功率谱密度的傅立叶逆变换
-
通过无偏信号的平方和将功率谱密度的逆傅里叶变换归一化,并且只取得一半结果vector
执行此操作的代码如下所示:
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)
我使用talib.CORREL进行自相关,我怀疑你可以用其他软件包做同样的事情:
def autocorrelate(x, period): # x is a deep indicator array # period of sample and slices of comparison # oldest data (period of input array) may be nan; remove it x = x[-np.count_nonzero(~np.isnan(x)):] # subtract mean to normalize indicator x -= np.mean(x) # isolate the recent sample to be autocorrelated sample = x[-period:] # create slices of indicator data correls = [] for n in range((len(x)-1), period, -1): alpha = period + n slices = (x[-alpha:])[:period] # compare each slice to the recent sample correls.append(ta.CORREL(slices, sample, period)[-1]) # fill in zeros for sample overlap period of recent correlations for n in range(period,0,-1): correls.append(0) # oldest data (autocorrelation period) will be nan; remove it correls = np.array(correls[-np.count_nonzero(~np.isnan(correls)):]) return correls # CORRELATION OF BEST FIT # the highest value correlation max_value = np.max(correls) # index of the best correlation max_index = np.argmax(correls)
我认为OP的问题的真正答案简洁地包含在Numpy.correlate文档的摘录中:
mode : {'valid', 'same', 'full'}, optional Refer to the `convolve` docstring. Note that the default is `valid`, unlike `convolve`, which uses `full`.
这意味着,当没有'模式'定义使用时,Numpy.correlate函数将返回一个标量,当给定两个input参数相同的向量(即 – 用于执行自相关)时。