周期检测算法详解：傅里叶变换、自相关系数与小波变换
周期检测算法详解：傅里叶变换、自相关系数与小波变换
周期检测是数据分析中的一个重要问题，可以帮助我们发现数据中的规律并预测未来的变化。本文将介绍三种常用的周期检测算法：傅里叶变换、自相关系数和小波变换。
背景
任何事物在两个不同时刻都不可能保持完全相同的状态，但很多变化往往存在着一定的规律，例如 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()
可以看到在单波、多波下，小波变换均有很好的效果