周期检测算法（自学使用，侵权立删）

背景

任何事物在两个不同时刻都不可能保持完全相同的状态，但很多变化往往存在着一定的规律，例如 24 小时日出日落，潮起潮落，这些现象通常称为「周期」。

周期性，指时间序列中呈现出来的围绕长期趋势的一种波浪形或振荡式变动。准确提取周期信息，不仅能反映当前数据的规律，应用于相关场景，还可以预测未来数据变化趋势。

时间序列示例

一般而言，时间序列周期性分为三种：

• 「符号性周期」，例如序列 fbcnfkgbfops中 f 周期为 4；

• 「部分周期性」，例如序列 ansdcdmncdcacdascdmc 中 cd 周期为 4；

• 「分段周期性」，例如上面给定的时间序列即为分段周期性；

针对时间序列的周期性检测问题，这篇文章主要介绍「傅里叶变换」和「自相关系数」两种方法及其在实际数据中的效果；

傅里叶变换

傅里叶变换是一种将时域、空域数据转化为频域数据的方法，任何波形（时域）都可以看做是不同振幅、不同相位正弦波的叠加（频域）。

傅里叶变换对于一条具备周期性的时间序列，它本身就很接近正弦波，所以它包含一个显著的正弦波，周期就是该正弦波的周期，而这个正弦波可以通过傅里叶变换找到，它将时序数据展开成三角函数的线性组合，得到每个展开项的系数，就是傅里叶系数。傅里叶系数越大，表明它所对应的正弦波的周期就越有可能是这份数据的周期。

自相关系数

自相关系数（Autocorrelation Function）度量的是同一事件不同时间的相关程度，不同相位差（lag）序列间的自相关系数可以用 Pearson 相关系数计算。当序列存在周期性时，遍历足够多的相位差，一定可以找到至少一个足够大的自相关系数，而它对应的相位差就是周期。所以对于检测时序周期来说，只需找到两个自相关系数达到一定阈值的子序列，它们起始时间的差值就是我们需要的周期。

为了保证结果的可靠性，可以将傅里叶分析和自相关系数结合起来判断周期性。主要思路是：先通过傅里叶变换找到可能的周期，再用自相关系数做排除，从而得到最可能的周期。

代码实例：

一、导入必要的库函数

import pywt
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False

 二、分别引入单波与合成波

sampling_rate = 1024 #采样率
t = np.arange(0, 1.0, 1.0 / sampling_rate) #采样区间

# 设定初始的四个频率
f1 = 20
f2 = 40
f3 = 50
f4 = 80
# 频率为20，振幅为1的单波函数

data = np.sin(2 * np.pi * f1 * t)

# 四个频率不同房、振幅相同的合成波函数
data = np.piecewise(t, [t < 1, t < 0.8, t < 0.5, t < 0.3],
                    [lambda t: 400*np.sin(2 * np.pi * f4 * t),
                     lambda t: 300*np.sin(2 * np.pi * f3 * t),
                     lambda t: 200*np.sin(2 * np.pi * f2 * t),
                     lambda t: 100*np.sin(2 * np.pi * f1 * t)])

三、 傅里叶变换

from scipy.fftpack import fft, fftfreq
 
fft_series = fft(data)
sample_freq = fftfreq(len(t), 1/sampling_rate)
fs = sample_freq[sample_freq>=0]
A = 2/sampling_rate*abs(fft_series[sample_freq>=0])
A[0] =  A[0]/2
 
plt.figure(figsize=(16,8))
plt.plot(fs,A,c='orangered',label='频率振幅')
plt.legend()
plt.grid(linestyle = ':')
plt.xlabel('Frequency(HZ)')
plt.ylabel('Amplitude(V)')
plt.title('傅里叶变换频率-幅值谱')

可以看出，在较为理想状况下，对于平稳的单波函数，傅里叶变换能够很好的把握周期，当然由于数据量不足，我们仅仅选用了正弦图像作为尝试，效果较好也可能与傅里叶变换本身的正弦、余弦映射相关联

四、 自相关系数验证

power = np.abs(fft_series)
pos_mask = np.where(sample_freq > 0)
powers = power[pos_mask]

top_k_seasons = 4
# 展示拟合程度最高的几个频率波
top_k_idxs = np.argpartition(powers, -top_k_seasons)[-top_k_seasons:]
top_k_power = powers[top_k_idxs]
fft_periods = (1 / fs[top_k_idxs]).astype(float)

print(f"top_k_power: {top_k_power}")
print(f"fft_periods: {fft_periods}")

from statsmodels.tsa.stattools import acf
fft_periods = (1 / fs[top_k_idxs]).astype(int)

# 期望整型的时间数据
for lag in fft_periods:
    acf_score = acf(data, nlags=lag)[-1]
    print(f"lag: {lag} fft acf: {acf_score}")

由于设定频率较高，T=0.2，不适宜用自相关系数图进行滞后判断，单已经可以看出结果与真实值的误差较小，仅为0.002*** 

傅里叶变换的局限性

      对于非平稳信号，傅里叶变换不能很好反映出其频率随时间的变化。因为我们在应用傅里叶变换时，计算出的每个频率分量都是对应于整个时间轴（或有信号的时间范围），这就使得原始信号的时间信息丢失了，不能分析出频率随时间的变化，也不能定位出某一时刻发生的突变。

      下图的三个信号在时域上有很大不同，第一张是不同频率信号相加混合在一起，后两张包含与第一张相同的频率成分，不过分时出现。但从傅里叶频谱上看，三者差别并不大，尤其是后两个非平稳信号，我们无法从频谱上区分它们。

 

 就拿我们上面的多波函数举例，FFT显然有点力不从心了

为了弥补FFT的不足，把整个时间域分解成无数个等长的小过程（加窗），每个过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率。

      但此时仍然有问题，那就是窗的宽度。窄窗口时间分辨率高、频率分辨率低；宽窗口时间分辨率低、频率分辨率高（如下图）。所以，对于时变非稳态信号，高频部分适合用小窗，低频部分适合用大窗。然而，在一次STFT中，窗口的宽度是固定。所以STFT也有其局限性，这就引出了我们的小波变换。

小波变换

 小波变换直接把傅里叶变换的基给换了，将无限长的三角函数基换成了有限长的会衰减的小波基，常用三角函数与指数函数的积作为小波基，它的能量有限，都集中在某一点附近，而且积分的值为零。傅里叶变换，变量只有w，而小波变换则有尺度a和平移量b，尺度对应于频率，平移量对应于时间，所以小波变换可以用于多个时段的时频分析，得到信号的时频谱。

我们选用了一个简单的小波基进行展示，并做了一定的可视化

wavename = 'cgau8'
totalscal = 256
fc = pywt.central_frequency(wavename)
cparam = 2 * fc * totalscal
scales = cparam / np.arange(totalscal, 1, -1)
[cwtmatr, frequencies] = pywt.cwt(data, scales, wavename, 1.0 / sampling_rate)

# 获取振幅信息
amplitude = np.abs(cwtmatr)

# 绘制原始信号
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(t, data)
plt.title('原始信号')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')

# 绘制振幅图
plt.subplot(212)
plt.imshow(amplitude, extent=[t.min(), t.max(), frequencies.min(), frequencies.max()], aspect='auto', cmap='PRGn', vmax=abs(amplitude).max(), vmin=-abs(amplitude).max())
plt.title('小波变换振幅图')
plt.xlabel('时间 (秒)')
plt.ylabel('频率 (Hz)')
plt.colorbar(label='振幅')

# 调整子图间距
plt.subplots_adjust(hspace=0.4)

# 显示图像
plt.show()

 可以看到在单波、多波下，小波变换均有很好的效果