Python信号处理：自相关函数(对标MATLAB中的autocorr)

摘要：Python中，更确切地说是numpy、scipy、statsmodels这些库中都有计算相关的方法。但numpy和scipy中的correlate方法的定义和MATLAB中的不同，导致计算结果不太一样。看上去MATLAB和statsmodels里都是用的标准的统计中的定义——皮尔森相关系数，而numpy和scipy中使用的是非正式的信号处理中的定义，需要均值为0，且计算结果需要归一化，才会得到差不多的答案。

系列目录

Python信号处理：快速傅里叶变换(FFT)，短时傅里叶变换(STFT)，窗函数，以及滤波
Python信号处理：自相关函数(对标MATLAB中的autocorr)
Python信号处理：波束形成及目标方位估计，CBF、MVDR
Python信号处理：cvxpy工具包求解稀疏约束优化问题

autocorr(series, 'NumLags', N-1)和statsmodels.tsa.api.stattools.acf(series, nlags=N-1)是相同的。[不考虑其他的参数]。
输入序列为0均值时，使用numpy.correlate(series, series, mode='full')，并对计算结果取后一半，得到的结果和前面两个方法基本相同。
输入序列不是0均值时，用numpy.correlate可能是不对的，因为信号处理下的非正式定义需要零均值条件，算出来的结果也很奇怪。

2. 非零均值序列的自相关函数

选择输入序列为y=[2, 3, 9, 6, 9]，N为5，分别使用上述三种方法得到自相关函数，得到的结果如下：
1.

acf = smt.stattools.acf(y, nlags=N-1)

array([ 1. , 0.06915888, -0.05794393, -0.2271028 , -0.28411215])

acf = np.correlate(y, y, mode='full')
acf = acf[N-1:]
acf = acf / acf[0]

array([1. , 0.66824645, 0.55450237, 0.18483412, 0.08530806])

acf = autocorr(y,'NumLags', N-1);

acf = 1.0000, 0.0692, -0.0579, -0.2271, -0.2841

可以看出，smt.stattools.acf(y, nlags=N-1)和autocorr(y,'NumLags', N-1)的结果是一样的，但np.correlate(y, y, mode='full')的结果差距较大，且所有的相关系数都是正的。

3. 零均值序列的自相关函数

仍使用输入序列y=[2, 3, 9, 6, 9]，N为5，但对其进行0均值化处理，得到的结果如下：

y = y - y.mean()
acf = smt.stattools.acf(y, nlags=N-1)

array([ 1. , 0.06915888, -0.05794393, -0.2271028 , -0.28411215])

和上一节的2.不同，但变得和1.、3.相同

y = y - y.mean()
acf = np.correlate(y, y, mode='full')
acf = acf[N-1:]
acf = acf / acf[0]

array([ 1. , 0.06915888, -0.05794393, -0.2271028 , -0.28411215])

和之前的3.是一样的

y = y - mean(y);
acf = autocorr(y,'NumLags', N-1);

acf = 1.0000, 0.0692, -0.0579, -0.2271, -0.2841

可以看出，对于零均值化的序列，三种方法的结果都是一样的。

4. 消除时滞的影响

当序列很长时，如下图，三种方法计算的自相关函数都会随着时间慢慢减小。在smt.stattools.acf中可以通过unbiased参数进行修正。对于np.correlate，可以用序列长度进行归一化，同样能起到消除时滞的效果。根据smt.stattools.acf文档的描述，两种做法好像是一样的。当然MATLAB中肯定也是可以这么操作的。

对无噪声的sin序列，分别计算自相关函数，看看时滞的影响，以及消除时滞后的效果。

# 左上，np.correlate，未消除时滞
acf1 = np.correlate(x, x, mode='full')
acf1 = acf1[N-1:]
acf1 = acf1 / acf1[0]
ax[0, 0].plot(t, acf1)# 右上，np.correlate，消除时滞
acf1 = np.correlate(x, x, mode='full')
acf1 = acf1[N-1:]
acf1 = acf1 / np.arange(N, 0, -1)
acf1 = acf1 / acf1[0]
ax[0, 1].plot(t, acf1)# 左下，smt.stattools.acf，未消除时滞
acf2 = smt.stattools.acf(x, nlags=N-1, unbiased=False)
ax[1, 0].plot(t, acf2)# 右下，smt.stattools.acf，消除时滞
acf2 = smt.stattools.acf(x, nlags=N-1, unbiased=True)
ax[1, 1].plot(t, acf2)