背景
任何事物在两个不同时刻都不可能保持完全相同的状态,但很多变化往往存在着一定的规律,例如 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()
可以看到在单波、多波下,小波变换均有很好的效果