计算时序数据的周期性

构造时序数据

import numpy as np
import matplotlib.pyplot as plt

# 设置参数
period = 128
num_cycles = 5
total_length = period * num_cycles

# 生成周期性信号（正弦波形）
np.random.seed(42)
time = np.arange(0, total_length, 1)
signal = 10 * np.sin(2 * np.pi * time / period)

# 加上噪声
noise = np.random.normal(0, 1, len(time))
signal_with_noise = signal + noise

# 打印生成的数据列表
signal = signal_with_noise.tolist()
# print(signal)

# 绘制周期性信号
plt.figure(figsize=(14, 4))
plt.plot(time, signal_with_noise, label='Sinusoidal Signal with Noise')
plt.xlabel('Time')
plt.ylabel('Value')
plt.title('Generated Periodic Signal with Noise')
plt.legend()
plt.show()

方法一：自相关函数

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
from scipy.signal import find_peaks

# 计算自相关函数
nlags = len(signal)/2  # 数据长度的一半
acf_values = acf(signal, nlags=nlags)

# 寻找显著峰值
peaks, _ = find_peaks(acf_values, height=0.5)  # height 参数决定阈值，可以调整

# 拟合周期（根据发现的峰值)
if len(peaks) > 1:
    peak_intervals = np.diff(peaks)  # 计算相邻峰值之间的距离
    estimated_period = np.mean(peak_intervals)
    print(f"Estimated period: {estimated_period:.2f} data points")
else:
    print("No significant period found. Consider adjusting the peak detection threshold.")

# 绘制自相关函数图
plt.figure(figsize=(14, 4))
plt.stem(acf_values, use_line_collection=True, markerfmt='C0o')
plt.plot(peaks, acf_values[peaks], "x", label='Significant Peaks', color='red')
plt.xlabel('Lags')
plt.ylabel('ACF')
plt.title('Autocorrelation Function with Detected Peaks')
plt.legend()
plt.show()

方法二：傅里叶变换

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftfreq


# 计算傅里叶变换
N = len(signal)
signal_fft = fft(signal)
frequencies = fftfreq(N, d=1)

# 仅保留非负频率部分
frequencies = frequencies[:N//2]
signal_fft = signal_fft[:N//2]

# 绘制频谱
plt.figure(figsize=(14, 4))
plt.stem(frequencies, np.abs(signal_fft), use_line_collection=True)
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.title('Frequency Domain (FFT)')
plt.tight_layout()
plt.show()

# 找到频率域中的峰值
peak_freq = frequencies[np.argmax(np.abs(signal_fft))]
estimated_period = 1 / peak_freq
print(f"Estimated period: {estimated_period:.2f} data points")