从时序数据中提取特征分量

 原始数据

import matplotlib.pyplot as plt
from matplotlib import font_manager
fname="/usr/local/python3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf"
zhfont = font_manager.FontProperties(fname=fname)


plt.figure(figsize=(14, 4))
plt.plot(dataset)
plt.title("原始数据", fontproperties=zhfont, fontsize=12)
plt.grid()
plt.show()

方法一:傅里叶变换

import datetime
import numpy as np
import pandas as pd
from scipy.fft import fft, ifft
import matplotlib.pyplot as plt
from matplotlib import font_manager
fname="/usr/local/python3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf"
zhfont = font_manager.FontProperties(fname=fname)


# 样本:288个5分钟点,即一天
N = len(dataset)

# 采样间隔(秒)
T = 300  # 每5分钟一个数据点

# 采样频率(Hz)
Fs = 1 / T  # 即每300秒一个样本

# 总观察时间(秒)
D = N * T  # 一整天的时间长度

# 生成时间向量
t = np.arange(0, D, T)

# 生成示例时序数据(两个正弦波的叠加)
# 假设每天有一个主要周期为24小时(低频)和一个次要周期为12小时(高频)
# y = 1.0 * np.sin(2.0 * np.pi * (1 / (24*3600)) * t) + 0.5 * np.sin(2.0 * np.pi * (1 / (12*3600)) * t)
y = dataset

# 计算傅里叶变换
yf = fft(y)
# 频率向量
xf = np.fft.fftfreq(N, T)

# 绘制频谱图
plt.figure(figsize=(14, 4))
plt.plot(xf, np.abs(yf))
plt.title("将原始数据从时域转成频域", fontproperties=zhfont, fontsize=12)
plt.xlabel("Frequency (Hz)", fontsize=12)
# plt.ylabel("Amplitude")
plt.grid()
plt.show()

# 提取各个正弦分量
def extract_sine_component(yf, freq_index):
    yf_component = np.zeros_like(yf)
    yf_component[freq_index] = yf[freq_index]
    if freq_index != 0:
        yf_component[-freq_index] = yf[-freq_index]  # 对称频谱
    y_component = ifft(yf_component)
    return np.real(y_component)

# 找出幅值最大的10个频率分量及其索引(只处理前N//2部分)
amplitudes = np.abs(yf[:N//2])
indices = np.argsort(amplitudes)[-10:]

plt.figure(figsize=(14, 4))
y_reconstructed = np.zeros_like(y)
for idx in indices:
    y_component = extract_sine_component(yf, idx)
    y_reconstructed += y_component
    plt.plot(t, y_component, label=f'Frequency: {xf[idx]:.5f} Hz')
plt.title("基于傅里叶变换将原始数据拆分成多个正弦波,取Top10", fontproperties=zhfont, fontsize=12)
plt.legend()
plt.grid()
plt.show()

# 绘制重建后的时序数据与原始数据的对比
plt.figure(figsize=(14, 4))
plt.plot(t, y, label='原始数据')
plt.plot(t, y_reconstructed, label='Top10正弦波叠加', linestyle='--')
plt.title("原始数据 vs Top10正弦波叠加", fontproperties=zhfont, fontsize=12)
plt.legend()
plt.grid()
plt.show()

 

posted @ 2024-07-04 18:27  lixin[at]hitwh  阅读(3)  评论(0编辑  收藏  举报