python 自相关序列(ACF)

学习目标（ILOS）：
您应该:
进一步熟悉相关性和自相关序列。
能够使用自相关序列来估计信号的基本自由度 /音高频率
能够产生震撼力
了解Python库的基本用途

# Let's do the ususal necessary and nice-to-have imports
%matplotlib notebook
import matplotlib.pyplot as plt     # plotting
import seaborn as sns; sns.set()    # styling ((un-)comment if you want)
import numpy as np                  # math

# imports we need in addition for this lab sheet
from IPython import display as ipd  # to listen to audio
import scipy.signal as sig          # for PSD calculation
import soundfile as sf              # to load WAVE files

L = 8000 # total length of signal 
P = 100  # period

K = 250  # analysis length

# generate periodic random signal
np.random.seed(1) # this ensures that always the same "random" signal is generated
x0 = np.random.normal(size=P)
x  = np.tile(x0, L//P)

plt.figure(figsize=(10, 4))
#plt.stem(x[:2*K], basefmt='C0:', use_line_collection=True)
plt.plot(x[:2*K])
plt.xlim(0, 2*K)
plt.xlabel('$k$')
plt.ylabel('$x[k]$')
plt.title('First '+str(2*K)+'samples of a periodically repeated noise signal')
plt.grid(True)

# listen to the sound file
ipd.Audio(x, rate=8000)

可以不设置参数，seed后random.random 返回的是一个任意的数字；如果设置参数后，只要参数不变，反复调用random.random方法（每调用一次该方法最好先运行random.seed()来产生新的随机数种子），他也只会返回一个相同的数字.
np.random.normal()的意思是一个正态分布，normal这里是正态的意思。参数loc(float)：正态分布的均值，对应着这个分布的中心。loc=0说明这一个以Y轴为对称轴的正态分布，
参数scale(float)：正态分布的标准差，对应分布的宽度，scale越大，正态分布的曲线越矮胖，scale越小，曲线越高瘦。
参数size(int 或者整数元组)：输出的值赋在shape里，默认为None。
棉棒图stem()x：制定棉棒的x轴基线上的位置y：绘制棉棒的长度
linefmt：棉棒的样式;markerfmt：棉棒末端的样式;basefmt：指定基线的样式

自相关序列，又称自相关功能
在上一个中，我们计算了相关系数和相关系数r^xx作为标量数量，指示相关的两个信号（相同长度）是如何相关的。 通常，相关性是对随机过程之间或一个随机过程样本之间依赖关系的统计度量。 结果，我们对X [K]和X [K-κ]之间的（自动）相关性感兴趣，X [K]和X [K]的时移版本感兴趣。 可以通过自动相关函数（ACF）来分析此信息，该函数表征了一个随机信号x [k]中的时间依赖性。 我们将使用ACF来估计本语音激发信号的基本（音高）频率。

背景理论和定义
如果我们假设广义平稳性（WSS），ACF确实取决于时间差κ= k1-k2在所考虑的样品索引和信号的标准偏差之间为σx=σy= 1，我们可以简化
rxx[κ]=E{x[k]⋅x[k−κ]} 其中通常选择κ作为样品指数而不是k表示它表示时间移位 /滞后。 ACF量化了信号与自身移动版本的相似性。 它具有高（正 /负）的高相似性和低相似性（正 /负）值的高（正 /负）值，这与我们的发现一致

相关估计器
由于取决于期望运算符E{⋅}，这只能估计而不是计算，我们对这样的估计器感兴趣 RXX [κ]。 因此，我们定义了在时间移位的时间移位κ的x [k]长度为L的时间移位

请注意，在数值实现（例如MATLAB，具有模式='full’的python）中，计算的ACF返回一个长度为2L -1的向量，其中包含负值从负κs到正κs。 对于关于实际时间滞后κ的解释，必须从这些指数中减去L -1。 更明显地，对于实值信号，自相关是对称的κ= 0。

进一步指出，eq：autocortrolationestfullsum 中的1L平均会产生ACF的偏置估计量，该估计量应始终用r^xx表示，偏置[κ]。 偏置是三角形的形式，可以通过适应归一化来去除。 ACF的公正估计器被给予

K = 250  # upper/lower limit for lag in ACF

# compute the ACF
acf = 1/len(x) * np.correlate(x, x, mode='full')

print('Length of the signal x is L=' + str(len(x)))
print('Length of the ACF is : ' + str(len(acf)))
print(acf)

# plot ACF over full length
kappa = np.arange(-(L-1), L)
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(kappa,acf)
plt.xlabel('time shift $\kappa$')
plt.ylabel('$\hat{r}_{xx}[\kappa]$')

# truncate ACF and plot area of interest around the center
acf_center=(len(x)-1)
kappa = np.arange(-(L-1), L)
plt.subplot(1,2,2)
plt.plot(kappa[acf_center-(K-1):acf_center+(K-1)],
         acf[acf_center-(K-1):acf_center+(K-1)])
plt.xlabel('time shift $\kappa$')
plt.ylabel('$\hat{r}_{xx}[\kappa]$')
None # suppress last output

np.correlate函数:numpy.correlate(a, v, mode=‘valid’)计算序列a,v的互相关。

我们可视化语音信号S [k]的开头，这是以下持续的vovel。 我们看到了周期性细分市场的开始。 以下图的下面板计算ACF r^ss [κ]。

K=1000

# compute and truncate ACF
acf = 1/len(s) * np.correlate(s, s, mode='full')
acf = acf[(len(s)-1)-(K-1):(len(s)-1)+K]
kappa = np.arange(-(K-1), K)
print(str(np.min(kappa)) + ' <= kappa <= ' + str(np.max(kappa)) )

# plot signal and its ACF
plt.figure(figsize=(10, 4))
#plt.stem(x[:4*K], basefmt='C0:', use_line_collection=True)
plt.plot(s[:4*K])
plt.xlim(0, 4*K)
plt.xlabel('$k$')
plt.ylabel('$s[k]$ for the first $2K='+str(4*K)+'$ samples')
plt.grid(True)

plt.figure(figsize=(10, 4))
#plt.stem(kappa, acf, basefmt='C0:', use_line_collection=True)
plt.plot(kappa, acf)
plt.xlabel('$\kappa$')
plt.ylabel('$\hat{r}_{ss}[\kappa]$  for ' + str(np.min(kappa)) + ' $\leq \kappa \leq ' + str(np.max(kappa)) + '$')
plt.axis([-K, K, 1.1*min(acf), 1.1*max(acf)])
plt.grid(True)

基本频率估计
借助ACF，我们现在将提取语音信号的基本频率F0，即音调频率，是语音生产系统的频率迁移，即声带之后。

为了确定F0，我们首先计算信号的自相关估计值R^XX [κ]，然后寻找最高的相关值。 为了确定信号的周期性，我们对ACF最大值的距离感兴趣。 由于我们知道我们始终在κ= 0时具有最大值，因此我们有兴趣识别第二个主要最大值。

在混合（发声/未发音的）激发信号的情况下更加健壮，我们只能在时间移位范围内看一下，这与人类声音的可能频率相对应。 人类的典型音高频率在FMIN = 50 Hz和FMAX = 500 Hz之间。 鉴于我们可以通过采样频率fs fs离散时间变量κ与离散频率变量有关的知识，我们可以在κMax= round（fs/fmax）和κmin= round（fs/fmin）之间计算搜索区域。

在这两个值之间，我们发现最大值的时变κ0为我们提供了相应的频率F0 =FSκ0。
该关系还显示在下图中，我们看到一些自相关序列，并观察FMIN和FMAX如何限制搜索空间的最高自相关值。

f_min_hz=50    # lower frequency in Hz defining our search range
f_max_hz=500   # upper frequency in Hz defining our search range

# The minimal and maximal shift (in samples) we want to look at is calculated from the frequency boundaries
kappa_acf_center=len(acf)//2
kappa_min = kappa_acf_center + int(np.round(fs / f_max_hz))
kappa_max = kappa_acf_center + int(np.round(fs / f_min_hz))

print('Length of ACF:    '+str(len(acf)))
print('ACF center index: '+str(kappa_acf_center))
print('kappa_min:        '+str(kappa_min))
print('kappa_max:        '+str(kappa_max))

# find the maximum value of the ACF (in the search range)
max_correlation_kappa = kappa_min + np.argmax(acf[kappa_min : kappa_max + 1])

# calculate pitch frequency
f_p=fs/(max_correlation_kappa-kappa_acf_center)

print('maximum at sample '+str(max_correlation_kappa)+', i.e. '+
      str(max_correlation_kappa-kappa_acf_center)+' samples from \kappa=0.')
print('pitch frequency f_0 is '+str(f_p)+' Hz')


plt.figure(figsize=(10, 4))
#plt.stem(kappa, acf, basefmt='C0:', use_line_collection=True)
plt.plot(kappa, acf)
plt.xlabel('$\kappa$')
plt.ylabel('$\hat{r}_{xx}[\kappa]$  for ' + str(np.min(kappa)) + ' $\leq \kappa \leq ' + str(np.max(kappa)) + '$')
plt.axis([-K, K, 1.1*min(acf), 1.1*max(acf)])

# add a shaded area illustrating the "search range"
plt.axvspan(kappa_min-kappa_acf_center,kappa_max-kappa_acf_center, alpha=0.5, color='orange')

# add vertical lines to the plot to visualise found values
plt.axvline(kappa_min-kappa_acf_center, color='r', ls='--')
plt.axvline(kappa_max-kappa_acf_center, color='b', ls='--')
plt.axvline(max_correlation_kappa-kappa_acf_center, color='g', ls='--')

# add some annotation
plt.text(kappa_min-kappa_acf_center,np.max(acf),
         '$f_{\mathrm{max}}='+str(f_max_hz)+'$ Hz')
plt.text(kappa_max-kappa_acf_center,np.max(acf),
         '$f_{\mathrm{min}}='+str(f_min_hz)+'$ Hz')
plt.text(kappa_min-kappa_acf_center,np.max(acf)-0.001,
         '$\kappa_{\mathrm{min}}='+str(kappa_min-kappa_acf_center)+'$')
plt.text(kappa_max-kappa_acf_center,np.max(acf)-0.001,
         '$\kappa_{\mathrm{max}}='+str(kappa_max-kappa_acf_center)+'$')
plt.text(max_correlation_kappa-kappa_acf_center,np.max(acf),
         '$f_0$ =%.2f Hz' % f_p)
plt.grid(True)

axvspan函数功能：在坐标系添加垂直区域（矩形）
在axvspan函数中，xmin，xmax是必输参数，无默认值，不写会报错
传入xmin，xmax设置区域的水平方向上的起始值与终止值，可以看到区域在水平方向上的4-6范围内，而竖直的y轴方向上则贯穿全部。
设置y轴方向的最大、最小值，修改默认整个y轴的区域设置，指定在y轴的区间，ymin=0.25,ymax=0.5,则y轴方向下边界为：y轴方向最小值+y轴区间长度* ymin,即− 1 + ( 1 − ( − 1 ) ) ∗ 0.25 = − 0.5 -1+(1-(-1))*0.25=-0.5−1+(1−(−1))∗0.25=−0.5,同理可得y轴方向上边界值。
修改参考与区域的颜色透明度,颜色通过color参数指定，透明度通过alpha参数指定，alpha取值位于[0,1]之间，取值为0，区域完全透明，不可见；取（0，1]之间的值即可。
axvline:绘制垂直于x轴的水平参考线