信号的相关分析（自相关和互相关） | python代码实战

欢迎关注公众号《故障诊断与python学习》
代码位置：https://github.com/HappyBoy-cmd/fault_diagnosis_signal_processing
参考资料：
书籍：机械故障诊断及典型案例解析（第2版，时献江）
Python信号处理：自相关函数(对标MATLAB中的autocorr)

一、自相关函数

(1) 自相关函数的定义

已知时间函数$x(t)$，其自相关函数$R_x(\tau )$的定义为：

$R_x(\tau )={\lim }_{T\to \infty }\frac{1}{T}{\int }_0^{T}x(t)x(t+\tau )\mathrm{d}t$

$R_x(\tau )$主要用来描述$x(t)$与其自身延时 $\tau$时刻之后的$x(t+\tau )$相似程度，相似程度越高，相关值越大。$R_x(\tau )$的计算原理如图3.9 所示。

图3.9 自相关函数计算示意图

(2) 自相关函数的性质

① $R_x(\tau )$为实偶函数，即$R_x(\tau )=R_x(-\tau )$ 。

② 当$\tau =0$时， $R_x(\tau )$的值最大，即 $R_x(0)⩾R_x(\tau )$， 并等于信号的均方值 ${\Psi }_x^2$。
$R_x(0)=\lim \frac{1}{T}{\int }_0^{T}x(t)x(t+0)\mathrm{d}t=\lim \frac{1}{T}{\int }_0^T{x}^2(t)\mathrm{d}t={\Psi }_x^2$

③ 当$\tau \to \infty$时，可认为$x(t)$和$x(t+\tau )$之间无关，即：
$R_x(\tau \to \infty )\to {\mu }_x^2$
若$x(t)$的均值为0, 则$R_x(\tau \to \infty )\to 0$。
性质② 和③ 的图像如图3.10 所示。

图3.10 自相关函数图

④ 周期函数的自相关函数仍是同频率的周期函数。假设正弦函数$x(t)=x_0\sin (\omega t+\phi )$的初始相位$\phi$是一个随机变量，根据自相关函数的定义，可求的自相关函数为：

$R_x(\tau )=\frac{x_0^2}{2}\cos \omega \tau$

可见正弦函数的自相关函数是一个余弦函数，如图3.11所示。在$\tau =0$时具有最大值式$x_0^2/2$它保留了变量$x(t)$的幅值信息$x_0$和频率$\omega$信息，但丢掉了初始相位$\phi$信息。

图3.11 正弦函数及其自相关函数

自相关函数主要用于检测混淆在随机信号中的周期信号成分，因为周期信号的自相关函数会按原频率重复出现，而随机信号在时间位移$\tau$稍大时，由于自身的相乘消除作用，自相关函数很快趋于零（ 假设均值${\mu }_x=0$ ) ，如图3.12 所示。

设测得信号$y(t)=x(t)+n(t)$，其中$x(t)$为正弦信号，$n(t)$为噪声，$x(t)$与$n(t)$相互独立，则有：

$R_y(\tau )=R_x(\tau )+R_n(\tau )$

当$\tau$很大时， $R_n(\tau )\to 0$ ， 因此：
${\lim }_{\tau \to \infty }\frac{R_y(\tau )}{R_x(\tau )}=1$

即当$\tau$大到一定程度时，$y(t)$与$x(t)$完全相关。

实际工程信号多数是随机噪声和确定性周期信号的混合体，如图3.13 所示。一般情况下．周期信号和故障特征有关，随机噪声对诊断无用。此时，可以利用随机信号的自相关函数迅速衰减，而周期函数不衰减的特性。在自相关图的右侧部分测取信号的周期，也就是说，自相关函数是从干扰噪声中找出周期信号或瞬时信号的重要手段． 延长变量$\tau$的取值，就可将信号中的周期分量${\tau }_0$暴露出来。

图3.13 随机噪声加周期信号的自相关函数

(3) 自相关函数的计算

由于采样点数的限制，按照式（3.31）进行自相关函数的计算是不可能的，因为我们只能得到有限的样本曲线及有限的数据长度。对于连续的模拟信号$x(t)$，如测量时间长度为 , 则其自相关函数可按下式计算：

${\hat{R}}_x(\tau )=\frac{1}{T-\tau }{\int }_0^{T-\tau }x(t)x(t+\tau )\mathrm{d}t(0⩽t⩽T)$

式中，${\hat{R}}_x(\tau )$ 表示$R_x(\tau )$的估计值。 时延$\tau$一定要远小于$T$, 以保证测量精度。

对于从连续信号采样所得的离散的数字信号$x(n),n=1,2,\cdots ,N$，其自相关函数可按下式估算：
${\hat{R}}_x(k)=\frac{1}{N-k}{\sum }_{n=1}^{N-k}x(n)(n+k)k=0,1,2,\cdots ,M;M\ll N$

为了保证测量精度，同样要使最大计算时延量$M$远远小于数据点数$N$, 以上可由计算机实现，称为直接计算法。因其计算量很大。近代的信号分析中已不采用这种方法。而是利用自相关函数与功率谱密度函数的关系，采用快速傅里叶变换算法实现，这部分内容将在后续介绍。

图 3.14 自相关函数计算结果图

二、互相关函数

(1) 互相关函数的定义

已知两个不同的信号$x(t)$和$y(t)$，$x(t)$和 $y(t)$的互相关函数$R_{xy}(\tau )$定义为：

$R_{xy}(\tau )={\lim }_{T\to \infty }\frac{1}{T}{\int }_0^{T}x(t)y(t+\tau )\mathrm{d}t$

互相关函数用于评价两个信号之间的相似程度，其计算原理如图3.15 所示。

图3.15 互相关函数计算示意图

(2) 互相关函数的性质

① 互相关函数是非奇、非偶函数， 即 $R_{xy}(\tau )=R_{yx}(-\tau )$。

② $\left|R_{xy}(\tau )\right|⩽\sqrt{R_{xx}(0)R_{yy}(0)}$，即$R_{xy}(0)$ 一般不是最大值，$R_{xy}(\tau )$的峰值不在$\tau =0$处。$R_{xy}(\tau )$的峰值偏离原点的位置${\tau }_0$反映了两信号时移的大小，此时两信号的相关程度最高，如图3.16 所示。
例如，当$x(t)=x_0\sin (\omega t+\theta )$，$y(t)=y0\sin (\omega t+\theta +\phi )$ 时， 其互相关函数$R_{xy}(\tau )$为：

$\begin{array}{l}{R}_{xy}(\tau )={\lim }_{T\to \infty }\frac{1}{T}{\int }_0^{T}x(t)y(t+\tau )\mathrm{d}t\\ =\frac{1}{T_0}{\int }_0^{T_0}{x}_0{y}_0\sin (\omega t+\theta )\sin [\omega (t+\tau )+\theta -\phi ]\mathrm{d}t\\ =\frac{x_0{y}_0}{2}\cos (\omega \tau -\phi )\end{array}$

由此可见，在$\tau =\phi /\omega ={\tau }_0$处，$R_{xy}(\tau )$取得最大值，表明两个信号此时最相关，因为$y(t)$延时$y(t)$时刻后与$x(t)$最相近。与自相关函数不同，两个同频率的谐波信号的互相关函数不仅保留了两个信号的幅值$x_0$，$y_0$信息、频率 $\omega$信息，而且还保留了两信号的相位差（ 此处为$\phi$）信息。
③ 两个统计独立的随机信号， 当均值为零时，则 $R_{xy}(\tau )=0$。
④ 两个不同频率的周期信号互不相关，即互相关函数为零。假设对两个不同频率的简谐波信号：
$x(t)=A_0\sin \left({\omega }_{1}t+\theta \right),y(t)=B_0\sin \left({\omega }_{2}t+\theta +\phi \right)$
进行相关分析，则：
R x y ( τ ) = lim ⁡ T → ∞ 1 T ∫ 0 T 0 x ( t ) y ( t + τ ) d t = 1 T 0 ∫ 0 T 0 A 0 B 0 sin ⁡ ( ω 1 t + θ ) sin ⁡ [ ω 2 ( t + τ ) + θ − φ ] d t = A 0 B 0 2 T 0 ∫ 0 T 0 { cos ⁡ [ ( ω 2 − ω 1 ) t + ( ω 2 τ − φ ) ] − cos ⁡ [ ( ω 2 + ω 1 ) t + ( ω 2 τ + 2 θ − φ ) ] } d t = 0 \begin{array}{l} R_{x y}(\tau)=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T_{0}} x(t) y(t+\tau) \mathrm{d} t \\ =\frac{1}{T_{0}} \int_{0}^{T_{0}} A_{0} B_{0} \sin \left(\omega_{1} t+\theta\right) \sin \left[\omega_{2}(t+\tau)+\theta-\varphi\right] \mathrm{d} t \\ =\frac{A_{0} B_{0}}{2 T_{0}} \int_{0}^{T_{0}}\left\{\cos \left[\left(\omega_{2}-\omega_{1}\right) t+\left(\omega_{2} \tau-\varphi\right)\right]-\cos \left[\left(\omega_{2}+\omega_{1}\right) t+\left(\omega_{2} \tau+2 \theta-\varphi\right)\right]\right\} \mathrm{d} t \\ =0 \end{array} Rxy​(τ)=limT→∞​T1​∫0T0​​x(t)y(t+τ)dt=T0​1​∫0T0​​A0​B0​sin(ω1​t+θ)sin[ω2​(t+τ)+θ−φ]dt=2T0​A0​B0​​∫0T0​​{cos[(ω2​−ω1​)t+(ω2​τ−φ)]−cos[(ω2​+ω1​)t+(ω2​τ+2θ−φ)]}dt=0​
即$R_{xy}(\tau )=0$ 。
⑤ 周期信号与随机信号的互相关函数为零。由于随机信号$y(t+\tau )$在时间$t\to t+\tau$内并无确定的关系，它的取值显然与任何周期函数$x(t)$无关，因此$R_{xy}(\tau )=0$ 。

三、相关分析的应用

(1) 自相关函数的应用

用轮廓仪测得一机械加工表面的粗糙度信号$a(t)$，并进行自相关分析，得到自相关函数$R_a(\tau )$，如图3.17 所示。根据$R_a(\tau )$就可以判断造成机械加工表面的粗糙度的原因。

图3.17 表面粗糙度的相关检测法

观察$a(t)$的自相关函数$R_a(\tau )$，发现其呈周期性，这说明造成粗糙度的原因之一是某种周期因素。从自相关函数图可以确定周期因素的频率为：
$f=\frac{1}{T}=\frac{1}{0.5/3}=6(\mathrm{Hz})$

根据加工该工件的机械设备中的各个运动部件的运动频率（如电动机的转速，拖板的往复运动次数，液压系统的油脉动频率等），通过测算和对比分析，运动频率与6Hz 接近的部件的振动，就应该是造成该粗糙度的主要原因。

(2) 互相关函数的应用

相关分析在机械设备故障诊断和振动控制中最直接的应用是传递问题，其中包括传递路径的识别和故障源的识别这两类问题。 间接的应用是相关测速和相关定位问题。

设时间信号（振动或噪声）$x(t)$通过一个非频变线性路径进行传递，传递过程中产生时延 ${\tau }_1$，并混入噪声$n(t)$。 可以用图3.18 (a) 描述。若传递路径的衰减因子为常数$\alpha$，则这个系统的输出$y(t)$可表示为：
$y(t)=\alpha x\left(t-{\tau }_1\right)+n(t)$
计算输入与输出信号的互相关函数为：
$\begin{array}{l}{R}_{xy}(\tau )={\lim }_{T\to \infty }\frac{1}{T}{\int }_0^{T}x(t)y(t+\tau )\mathrm{d}t\\ ={\lim }_{T\to \infty }\frac{1}{T}{\int }_0^{T}x(t)\left[\alpha x\left(t+\tau -{\tau }_1\right)+n(t+\tau )\right]\mathrm{d}t\\ ={\lim }_{T\to \infty }\frac{1}{T}{\int }_0^T\alpha x(t)x\left(t+\tau -{\tau }_1\right)+{\lim }_{T\to \infty }\frac{1}{T}{\int }_0^{T}x(t)n(t+\tau )\mathrm{d}t\\ =R_x\left(\tau -{\tau }_1\right)+R_{xn}(\tau )=R_x\left(\tau -{\tau }_1\right)\end{array}$、

根据互相关的性质⑤ ，式中$R_{xn}(\tau )=0$, 计算结果为$x(t)$在$\tau ={\tau }_1$的互相关函数，如图3.18（b）所示。该图表示$y(t)$在延时${\tau }_1$时刻后与$x(t)$相关．或者说$y(t)$是延时${\tau }_1$时刻后，并且幅值衰减了的$x(t)$。

图3.18 某信号的传递过程

① 故障源检测图3.19 给出一个用互相关函数诊断汽车驾驶员座椅振动源的例子。座椅上的振动信号为$y(t)$，前轮轴梁和后轮轴架上的振动信号分别为$x(t)$和$z(t)$，分别求$R_{xy}(\tau )$和$R_{zy}(\tau )$。从图上可看出$R_{xy}(\tau )$有突出的谱峰．说明座椅的振动$y(t)$主要是由于前轮的振动$x(t)$引起的。

图3.19 利用相关分析进行故障源检测

② 相关定位用
互相关分析法确定深埋地下的输油管裂损位置，如图3.20 所示。漏损处$K$可视为向两侧传播声音的声源，在两侧管道上分别放置传感器1 和2 。因为放置传感器的两点相距漏损处距离不等， 则漏油的声响传至两传感器的时间就会有差异，在互相关函数图上$\tau ={\tau }_{\mathrm{m}}$处有最大值，这个${\tau }_{\mathrm{m}}$就是时差，用${\tau }_{\mathrm{m}}$就可以确定漏损处的位置。

假设$v$为声音在管道中的传播速度，则漏损处$K$ 距中心点$O$ 的距离$S$ 为：
$S=\frac{1}{2}v{\tau }_{\mathrm{m}}$

图3.20 利用相关分析进行线性定位实例

四、python代码实战

4.1 自相关代码实战

4.1.1 导入包

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, fftpack, stats
import random
import os
import statsmodels.tsa.api as smt
import statsmodels

from matplotlib import rcParams

config = {
    "font.family": 'serif', # 衬线字体
    "font.size": 10, # 相当于小四大小
    "font.serif": ['SimSun'], # 宋体
    "mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    'axes.unicode_minus': False # 处理负号，即-号
}
rcParams.update(config)

fs = 1000   # 采样频率1000Hz
f = 10      # 信号频率 10Hz
n = np.arange(0, 4096)  # 生成4096点的序列
print(n)

>>>输出结果
[   0    1    2 ... 4093 4094 4095]

t = n/fs    # 生成4096点的时间序列
x = np.sin(2 * np.pi * n/fs * f) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize=(12,2))
plt.plot(x)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('正弦信号+噪声信号图')
plt.show()
print(x.shape)

>>>输出结果
(4096,)

4.1.2 绘制不同信号的自相关信号

plt.figure(figsize=(12,10))   
plt.subplots_adjust(wspace = 1, hspace = 1)
##===========绘制正弦信号图==========##
x1 = np.sin(2 * np.pi * n/fs * f)
plt.subplot(6, 1, 1)
plt.plot(t, x1)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('sin(2 * pi * t * f) 正弦信号图')

##===========绘制正弦信号+噪声图==========##
x2 = np.sin(2 * np.pi * n/fs * f) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.subplot(6, 1, 2)
plt.plot(t, x2)
plt.plot(t, x2)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('sin(2 * pi * t * f) 正弦信号+噪声图')

##===========绘制噪声信号图==========##
x3 = 0.5*np.random.randn(4096)
plt.subplot(6, 1, 3)
plt.plot(t, x3)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('噪声信号图')

##===========绘制正弦信号的自相关信号图==========##
acf1 = np.correlate(x1, x1, mode='full')
N = len(x1)
acf1 = acf1[N-1:]
acf1 = acf1 / np.arange(N, 0, -1)
acf1 = acf1 / acf1[0]
plt.subplot(6, 1, 4)
plt.plot(t, acf1)
plt.plot(t, acf1)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('正弦信号的自相关信号图')

##===========绘制正弦信号+噪声的自相关信号图==========##
acf2 = np.correlate(x2, x2, mode='full')
N = len(x2)
acf2 = acf2[N-1:]
acf2 = acf2 / np.arange(N, 0, -1)
acf2 = acf2 / acf2[0]
plt.subplot(6, 1, 5)
plt.plot(t, acf2)
plt.plot(t, acf2)
plt.plot(t, acf2)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('正弦信号+噪声的自相关信号图')


##===========绘制噪声信号的自相关信号图=========##
x4 = 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
acf3 = np.correlate(x4, x4, mode='full')
N = len(x4)
acf3 = acf3[N-1:]
acf3 = acf3 / np.arange(N, 0, -1)
acf3 = acf3 / acf3[0]
plt.subplot(6, 1, 6)
plt.plot(t, acf3)
plt.plot(t, acf3)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('噪声的自相关信号图')

由图可知，正弦信号的自相关信号为余弦函数。

4.1.3 自相关代码分析

x4 = np.sin(2 * np.pi * n/fs * f)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize=(12,6))   
plt.subplots_adjust(wspace = 1, hspace = 1)
##===========绘制正弦信号自相关信号未处理图=========##
plt.subplot(3, 1, 1)
acf3 = np.correlate(x4, x4, mode='full')
plt.plot(acf3)
##===========绘制正弦信号自相关信号一半图==========##
plt.subplot(3, 1, 2)
N = len(x4)
acf3 = acf3[N-1:]
plt.plot(acf3)
##===========绘制正弦信号的自相关信号处理后图=========##
plt.subplot(3, 1, 3)
acf3 = acf3 / np.arange(N, 0, -1)  # 用于消除时滞的影响 https://blog.csdn.net/qq_38290475/article/details/106992733
acf3 = acf3 / acf3[0]
plt.plot(acf3)

4.2互相关代码实战

4.2.1 同频率正弦余弦信号的互相关信号

t = n/fs    # 生成4096点的时间序列
f1 = 10
x1 = np.sin(2 * np.pi * n/fs * f1) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize = (12, 2))
plt.plot(x1)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('频率为10Hz正弦信号+噪声的时域图')

t = n/fs    # 生成4096点的时间序列
f2 = 10
x2 = 1*np.cos(2 * np.pi * n/fs * f2) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize = (12, 2))
plt.plot(x2)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('频率为10Hz余弦信号+噪声的时域图')

acf3 = np.correlate(x1, x2, mode='full')
N = len(x1)
acf3 = acf3[N-1:]
acf3 = acf3 / np.arange(N, 0, -1)
acf3 = acf3 / acf3[0]
plt.figure(figsize = (12, 2))
plt.plot(t, acf3)
plt.plot(t, acf3)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('同频率正弦余弦信号的自相关信号图（sin，cos）')

acf3 = np.correlate(x2, x1, mode='full')
N = len(x1)
acf3 = acf3[N-1:]
acf3 = acf3 / np.arange(N, 0, -1)
acf3 = acf3 / acf3[0]
plt.figure(figsize = (12, 2))
plt.plot(t, acf3)
plt.plot(t, acf3)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('同频率正弦余弦信号的自相关信号图（cos，sin）')

在进行自相关时，做自相关信号的前后顺序也会影响结果的。

4.2.2 不同频率信号之间的互相关信号

t = n/fs    # 生成4096点的时间序列
f1 = 10
x1 = np.sin(2 * np.pi * n/fs * f1) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize = (12, 2))
plt.plot(x1)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('频率为10Hz的正弦信号+噪声信号的时域图')

t = n/fs    # 生成4096点的时间序列
f2 = 15
x2 = np.cos(2 * np.pi * n/fs * f2) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize = (12, 2))
plt.plot(x2)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('频率为20Hz余弦信号+噪声信号的时域图')

acf3 = np.correlate(x1, x2, mode='full')
N = len(x1)
acf3 = acf3[N-1:]
acf3 = acf3 / np.arange(N, 0, -1)
acf3 = acf3 / acf3[0]
plt.figure(figsize = (12, 2))
plt.plot(t, acf3)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('不同频率正弦和余弦信号的自相关信号图')

4.2.3 正弦与噪声之间的互相关信号

t = n/fs    # 生成4096点的时间序列
f1 = 10
x1 = 0.5*np.random.randn(4096)  # 生成随机信号
plt.figure(figsize = (12, 2))
plt.plot(x1)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('噪声信号的时域图')
x.shape

t = n/fs    # 生成4096点的时间序列
f2 = 15
x2 = 1*np.cos(2 * np.pi * n/fs * f2) + 0.5*np.random.randn(4096)  # 生成正弦信号 sin(2*pi*f*t) + 随机信号
plt.figure(figsize = (12, 2))
plt.plot(x2)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('余弦信号+噪声的时域图')

acf3 = np.correlate(x1, x2, mode='full')
N = len(x1)
acf3 = acf3[N-1:]
acf3 = acf3 / np.arange(N, 0, -1)
acf3 = acf3 / acf3[0]
plt.figure(figsize = (12, 2))
plt.plot(t, acf3)
plt.plot(t, acf3)
plt.xlabel('time(s)')
plt.ylabel('Amp(m/s2)')
plt.title('正弦信号与噪声的互相关信号图')

五、结论

自相关

• 周期信号做自相关可以滤去噪声
• 噪声信号做自相关结果趋于0
• 未知周期信号+噪声信号做自相关可以获取周期信号的周期

互相关

• 同频率周期信号做互相关可以获取相位差
• 不同频率周期信号做互相关结果趋于0
• 周期信号与噪声信号做互相关结果趋于0