当前位置:网站首页>数字信号处理-8-自相关

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

2022-11-24 21:47:14小何才露尖尖角

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数字信号处理应用

原网站

版权声明
本文为[小何才露尖尖角]所创,转载请带上原文链接,感谢
https://blog.csdn.net/weixin_40994552/article/details/127996802

随机推荐