数字信号处理-8-自相关

1 皮尔森相关系数

假设 x 和 y 均为 N 个样本的数组，皮尔森公式如下：

皮尔森相关系数总是在 -1 到 +1 之间（包含这两个字）。ρ 的绝对值意味着相关性的强度。ρ 接近 +1 表示强正相关；ρ 接近 -1 表示强负相关，即随着一个值的增大另一个值减小。

如计算两个相位差为 1 的 sin 函数的相关性，从图形中可以看出两者具有相关性，一个升高，另一个也升高：

皮尔森相关系数矩阵如下，两者相关性为 0.52

皮尔森相关性计算代码：

# jupyter notebook 中运行
# 导入需要的包
import numpy as np
from matplotlib import pyplot as plt
import math

# 常数值 2π
# jupyter notebook 中运行
# 导入需要的包
import numpy as np
from matplotlib import pyplot as plt
import math

# 常数值 2π
PI2 = math.pi * 2
#采样率
framerate = 11025
# 采样时长：秒
duration = 0.01

def make_ys(amp=1, f=440, offset=0):
    """
    amp:振幅
    f:频率
    offset：相位偏移
    
    """
    # 样本数
    n = round(duration * framerate)
    # 采样时间点
    ts = np.arange(n) / framerate
    ys = amp * np.sin(PI2*f*ts + offset)
    return ys

wave1 = make_ys(offset=0)
wave2 = make_ys(offset=1)
# ddof=0 表示corrcoef 应该除以 N,而不是默认的 N-1
corr_matrix = np.corrcoef(ys1, ys2, ddof=0)
corr_matrix

2 序列相关

信号代表的是随时间变化的数值的度量。比如声音信号代表的是对电压的度量，对应于空气压力的变化，也即我们所感知到的声音。

这样的信号前后序列之间是相关的，为计算序列相关性，我们可以平移信号，然后计算平移之后的信号与原始信号之间的相关性。

# 上一节函数，产生序列
wave1 = make_ys(offset=0)
# 计算自相关
def serial_corr(wave, lag=1):
    # 序列长度
    n = len(wave)
    y1 = wave[lag:]
    y2 = wave[:n-lag]
    # 相关系数值
    corr = np.corrcoef(y1, y2, ddof=0)[0,1]
    return corr
serial_corr(wave1)

结果为 0.969655978695559

3 自相关

上一节只计算了 lag=1 的相关性，实际上可以设 lag 为不同值，从而得到一个不同 lag 对应的相关性值序列。

def autocorr(wave):
	# 最大 lag 为波形长度的一半，因为超过一半的话前后序列就不相交了
    lags = range(len(wave)//2)
    corrs = [serial_corr(wave, lag) for lag in lags]
    return lags, corrs
lags, corrs = autocorr(wave1)
plt.plot(lags, corrs)

lags 是一串从 0 到半个波形长的整数序列，corrs 是对应各个 lag 的序列相关性值的数组。

因为 sin 函数是周期函数，所以间隔一定 lag 之后，它的波形与原始波形又变为相同，相关性为 1。

据此在一段任意的波形中，通过计算序列自相关性找到最大的自相关位置，该位置我们就可以认为是此波形的周期。

3.1 通过自相关计算频率

如上图所示，自相关最大的第一个 lag 位置为 0（排除），随后 lag 滞后 25 相关性变最大，因此该点即为周期点。滞后 25 个样本之后波形变得重复，采样率为 11025，25个样本对应时间为 25/11025 秒，因此周期为 25/11025 秒。频率为周期的倒数，所以频率为 441Hz=11025/25。非常接近与我们设定的 440Hz。

3.2 Numpy 计算自相关

NumPy 提供了一个函数 correlate 来计算两个函数之间的相关性或者一个函数的自相关。

# mode='same' 对应的时滞范围为 -N/2 到 N/2,N为波形数组长度
corrs2 = np.correlate(wave1,wave1,mode='same')
N = len(corrs2)
x = range(-N//2, N//2, 1)
plt.plot(x, corrs2)

结果是对称的，这是因为针对同一个信号，一个信号的负时滞和对另一个结果的正时滞是一样的，所以取一半即可

# 时滞取一半，后一半
N = len(corrs2)
half = corrs2[N//2:]
# np.correlate 计算相关时，并没有除序列长度，可以通过除序列长度纠正
lengths = range(N, N//2, -1)
half /= lengths
# 结果归一化，这样 lag=0 时的相关性为 1
half /= half[0]
plt.plot(half)

与之前定义的函数结果相同 autocorr(wave)

参考

漫画傅里叶解析
Python数字信号处理应用