时间序列特征提取:从理论到Python代码实践

本文涉及的产品
智能开放搜索 OpenSearch行业算法版,1GB 20LCU 1个月
实时计算 Flink 版,5000CU*H 3个月
检索分析服务 Elasticsearch 版,2核4GB开发者规格 1个月
简介: 时间序列是一种特殊的存在。这意味着你对表格数据或图像进行的许多转换/操作/处理技术对于时间序列来说可能根本不起作用。

时间序列是一种特殊的存在。这意味着你对表格数据或图像进行的许多转换/操作/处理技术对于时间序列来说可能根本不起作用。

"特征提取"的想法是对我们拥有的数据进行"加工",确保我们提取所有有意义的特征,以便下一步(通常是机器学习应用)可以从中受益。也就是说它是一种通过提供重要特征并过滤掉所有不太重要的特征来"帮助"机器学习步骤的方法。

这是完整的特征提取过程:

对于表格数据和信号,他们的特征根本就不同,比如说的概念,傅里叶变换小波变换的想法,以及独立分量分析(ICA)的概念只有在处理信号时才真正有意义。

目前有两大类进行特征提取的方法:

  • 基于数据驱动的方法: 这些方法旨在仅通过观察信号来提取特征。它们忽略机器学习步骤及其目标(例如分类、预测或回归),只看信号,对其进行处理,并从中提取信息。
  • 基于模型的方法: 这些方法着眼于整个流程,旨在为特定问题找到解决方案的特征

数据驱动方法的优点是它们通常在计算上简单易用,不需要相应的目标输出。它们的缺点是特征不是针对你的特定问题的。例如,对信号进行傅里叶变换并将其用作特征可能不如在端到端模型中训练的特定特征那样优化。

在这篇文章中,为了简单起见,我们将只讨论数据驱动方法。并且讨论将基于领域特定的方法、基于频率的方法、基于时间的方法和基于统计的方法。

1、领域特定特征提取

提取特征的最佳方法是考虑特定的问题。例如,假设正在处理一个工程实验的信号,我们真正关心的是t = 6s后的振幅。这些是特征提取在一般情况下并不真正有意义,但对特定情况实际上非常相关。这就是我们所说的领域特定特征提取。这里没有太多的数学和编码,但这就是它应该是的样子,因为它极度依赖于你的特定情况。

2、基于频率的特征提取

这种方法与我们的时间序列/信号的谱分析有关。我们有一个自然域,自然域是查看信号的最简单方式,它是时间域,意味着我们将信号视为给定时间(或向量)。

例如考虑这个信号,在其自然域中:

绘制它就可以得到:

这是自然(时间)域,也是我们数据集的最简单域。

然后我们可以将其转换为频率域。 信号有三个周期性组件。频率域的想法是将信号分解为其周期性组件的频率振幅相位

信号y(t)的傅里叶变换Y(k)如下:

这描述了频率为k的分量的振幅相位。在提取有意义特征方面,可以提取10个主要分量(振幅最高的)的振幅、相位和频率值。这将是10x3个特征(振幅、频率和相位 x 10),它们将根据谱信息描述时间序列。

这种方法可以扩展。例如,可以将我们的信号分解,不是基于正弦/余弦函数,而是基于小波,这是另一种形式的周期波。这种分解被称为小波分解

我们使用代码来详细解释,首先从非常简单的傅里叶变换开始。

首先,我们需要邀请一些朋友来参加派对:

 importnumpyasnp
 importmatplotlib.pyplotasplt
 importpandasaspd

现在让我们以这个信号为例:

 plt.figure(figsize=(10,5))
 x=np.linspace(-2*np.pi,2*np.pi,1000)
 y=np.sin(x) +0.4*np.cos(2*x) +2*np.sin(3.2*x)
 plt.plot(x,y,color='firebrick',lw=2)
 plt.xlabel('Time (t)',fontsize=24)
 plt.xticks(fontsize=14)
 plt.yticks(fontsize=14)
 plt.grid(alpha=0.4)
 plt.ylabel('y',fontsize=24)

这个信号有三个主要组成部分。一个振幅 = 1和频率 = 1,一个振幅 = 0.4和频率 = 2,一个振幅 = 2和频率 = 3.2。我们可以通过运行傅里叶变换来恢复它们:

 importnumpyasnp
 importmatplotlib.pyplotasplt

 # Generate the time-domain signal
 x=np.linspace(-8*np.pi, 8*np.pi, 10000)
 y=np.sin(x) +0.4*np.cos(2*x) +2*np.sin(3.2*x)
 y=y-np.mean(y)
 # Perform the Fourier Transform
 Y=np.fft.fft(y)
 # Calculate the frequency bins
 frequencies=np.fft.fftfreq(len(x), d=(x[1] -x[0]) / (2*np.pi))
 # Normalize the amplitude of the FFT
 Y_abs=2*np.abs(Y) /len(x)
 # Zero out very small values to remove noise
 Y_abs[Y_abs<1e-6] =0
 relevant_frequencies=np.where((frequencies>0) & (frequencies<10))
 Y_phase=np.angle(Y)[relevant_frequencies]
 frequencies=frequencies[relevant_frequencies]
 Y_abs=Y_abs[relevant_frequencies]

 # Plot the magnitude of the Fourier Transform
 plt.figure(figsize=(10, 6))
 plt.plot(frequencies, Y_abs)
 plt.xlim(0, 10)  # Limit x-axis to focus on relevant frequencies
 plt.xticks([3.2,1,2])
 plt.title('Fourier Transform of the Signal')
 plt.xlabel('Frequency (Hz)')
 plt.ylabel('Magnitude')
 plt.grid(True)
 plt.show()

我们可以清楚地看到三个峰值,对应的振幅和频率。

可以用一个非常简单的函数来完成所有工作

 defextract_fft_features( y, x=None,  num_features=5,max_frequency=10):
   y=y-np.mean(y)
   # Perform the Fourier Transform
   Y=np.fft.fft(y)
   # Calculate the frequency bins
   ifxisNone:
     x=np.linspace(0,len(y))
   frequencies=np.fft.fftfreq(len(x), d=(x[1] -x[0]) / (2*np.pi))
   Y_abs=2*np.abs(Y) /len(x)
   Y_abs[Y_abs<1e-6] =0
   relevant_frequencies=np.where((frequencies>0) & (frequencies<max_frequency))
   Y_phase=np.angle(Y)[relevant_frequencies]
   frequencies=frequencies[relevant_frequencies]
   Y_abs=Y_abs[relevant_frequencies]
   largest_amplitudes=np.flip(np.argsort(Y_abs))[0:num_features]
   top_5_amplitude=Y_abs[largest_amplitudes]
   top_5_frequencies=frequencies[largest_amplitudes]
   top_5_phases=Y_phase[largest_amplitudes]
   fft_features=top_5_amplitude.tolist()+top_5_frequencies.tolist()+top_5_phases.tolist()
   amp_keys= ['Amplitude '+str(i) foriinrange(1,num_features+1)]
   freq_keys= ['Frequency '+str(i) foriinrange(1,num_features+1)]
   phase_keys= ['Phase '+str(i) foriinrange(1,num_features+1)]
   fft_keys=amp_keys+freq_keys+phase_keys
   fft_dict= {fft_keys[i]:fft_features[i] foriinrange(len(fft_keys))}
   fft_data=pd.DataFrame(fft_features).T
   fft_data.columns=fft_keys
   returnfft_dict, fft_data

所以如果输入信号y和(可选):

  • x时间数组,考虑的特征数量(或峰值),最大频率

输出:

 x=np.linspace(-8*np.pi, 8*np.pi, 10000)
 y=np.sin(x) +0.4*np.cos(2*x) +2*np.sin(3.2*x)
 extract_fft_features(x=x,y=y,num_features=3)[1]

如果想使用小波(不是正弦/余弦)*来提取特征,可以进行小波变换。需要安装这个库:

 pip install PyWavelets

然后运行这个:

 importnumpyasnp
 importpywt
 importpandasaspd

 defextract_wavelet_features(y, wavelet='db4', level=3, num_features=5):
     y=y-np.mean(y)  # Remove the mean

     # Perform the Discrete Wavelet Transform
     coeffs=pywt.wavedec(y, wavelet, level=level)

     # Flatten the list of coefficients into a single array
     coeffs_flat=np.hstack(coeffs)

     # Get the absolute values of the coefficients
     coeffs_abs=np.abs(coeffs_flat)

     # Find the indices of the largest coefficients
     largest_coeff_indices=np.flip(np.argsort(coeffs_abs))[0:num_features]

     # Extract the largest coefficients as features
     top_coeffs=coeffs_flat[largest_coeff_indices]

     # Generate feature names for the wavelet features
     feature_keys= ['Wavelet Coeff '+str(i+1) foriinrange(num_features)]

     # Create a dictionary for the features
     wavelet_dict= {feature_keys[i]: top_coeffs[i] foriinrange(num_features)}

     # Create a DataFrame for the features
     wavelet_data=pd.DataFrame(top_coeffs).T
     wavelet_data.columns=feature_keys

     returnwavelet_dict, wavelet_data

 # Example usage:
 wavelet_dict, wavelet_data=extract_wavelet_features(y)
 wavelet_data

3.、基于统计的特征提取

特征提取的另一种方法是依靠统计学。 给定一个信号,有多种方法可以从中提取一些统计信息。从简单到复杂,这是我们可以提取的信息列表:

  • 平均值,只是信号的总和除以信号的时间步数。

  • 方差,即信号与平均值的偏离程度:

  • 偏度峰度。这些是测试你的时间序列分布"不是高斯分布"程度的指标。偏度描述它有多不对称,峰度描述它有多"拖尾"。
  • 分位数:这些是将时间序列分成具有概率范围的区间的值。例如,0.25分位数值为10意味着你的时间序列中25%的值低于10,剩余75%大于10
  • 自相关。这基本上告诉你时间序列有多"模式化",意味着时间序列中模式的强度。或者这个指标表示时间序列值与其自身过去值的相关程度。
  • :表示时间序列的复杂性或不可预测性。

这些属性中的每一个都可以用一行代码实现:

 importnumpyasnp
 importpandasaspd
 fromscipy.statsimportskew, kurtosis

 defextract_statistical_features(y):
     defcalculate_entropy(y):
         # Ensure y is positive and normalized
         y=np.abs(y)
         y_sum=np.sum(y)

         # Avoid division by zero
         ify_sum==0:
             return0

         # Normalize the signal
         p=y/y_sum

         # Calculate entropy
         entropy_value=-np.sum(p*np.log2(p+1e-12))  # Add a small value to avoid log(0)

         returnentropy_value
     # Remove the mean to center the data
     y_centered=y-np.mean(y)
     y=y+np.max(np.abs(y))*10**-4

     # Statistical features
     mean_value=np.mean(y)
     variance_value=np.var(y)
     skewness_value=skew(y)
     kurtosis_value=kurtosis(y)
     autocorrelation_value=np.correlate(y_centered, y_centered, mode='full')[len(y) -1] /len(y)
     quantiles=np.percentile(y, [25, 50, 75])
     entropy_value=calculate_entropy(y)  # Add a small value to avoid log(0)

     # Create a dictionary of features
     statistical_dict= {
         'Mean': mean_value,
         'Variance': variance_value,
         'Skewness': skewness_value,
         'Kurtosis': kurtosis_value,
         'Autocorrelation': autocorrelation_value,
         'Quantile_25': quantiles[0],
         'Quantile_50': quantiles[1],
         'Quantile_75': quantiles[2],
         'Entropy': entropy_value
     }

     # Convert to DataFrame for easy visualization and manipulation
     statistical_data=pd.DataFrame([statistical_dict])

     returnstatistical_dict, statistical_data

 y=np.sin(x) +0.4*np.cos(2*x) +2*np.sin(3.2*x)
 wavelet_dict, wavelet_data=extract_statistical_features(y)
 wavelet_data

4、基于时间的特征提取器

这部分将专注于如何通过仅提取时间特征来提取时间序列的信息。比如提取峰值和谷值的信息。我们将使用来自SciPy的find_peaks函数

 importnumpyasnp
 fromscipy.signalimportfind_peaks

 defextract_peaks_and_valleys(y, N=10):
     # Find peaks and valleys
     peaks, _=find_peaks(y)
     valleys, _=find_peaks(-y)

     # Combine peaks and valleys
     all_extrema=np.concatenate((peaks, valleys))
     all_values=np.concatenate((y[peaks], -y[valleys]))

     # Sort by absolute amplitude (largest first)
     sorted_indices=np.argsort(-np.abs(all_values))
     sorted_extrema=all_extrema[sorted_indices]
     sorted_values=all_values[sorted_indices]

     # Select the top N extrema
     top_extrema=sorted_extrema[:N]
     top_values=sorted_values[:N]

     # Pad with zeros if fewer than N extrema are found
     iflen(top_extrema) <N:
         padding=10-len(top_extrema)
         top_extrema=np.pad(top_extrema, (0, padding), 'constant', constant_values=0)
         top_values=np.pad(top_values, (0, padding), 'constant', constant_values=0)

     # Prepare the features
     features= []
     foriinrange(N):
         features.append(top_values[i])
         features.append(top_extrema[i])

     # Create a dictionary of features
     feature_dict= {f'peak_{i+1}': features[2*i] foriinrange(N)}
     feature_dict.update({f'loc_{i+1}': features[2*i+1] foriinrange(N)})

     returnfeature_dict,pd.DataFrame([feature_dict])

 # Example usage:
 x=np.linspace(-2*np.pi,2*np.pi,1000)
 y=np.sin(x) +0.4*np.cos(2*x) +2*np.sin(3.2*x)
 features=extract_peaks_and_valleys(y, N=10)
 features[1]

我们选择N = 10个峰值,但信号实际上只有M = 4个峰值,剩余的6个位置和峰值振幅将为0。

5、使用哪种方法?

我们已经看到了4种不同类别的方法。

我们应该使用哪一种?

如果你有一个基于领域的特征提取,那总是最好的选择:如果实验的物理或问题的先验知识是清晰的,你应该依赖于它并将这些特征视为最重要的,甚至可能将它们视为唯一的特征。

频率统计基于时间的特征而言,在可以将它们一起使用。将这些特征添加到你的数据集中,然后看看其中的一些是否有帮助,没有帮助就把他们删除掉,这是一个测试的过程,就像超参数搜索一样,我们无法确定好坏,所以只能靠最后的结果来判断。

https://avoid.overfit.cn/post/5790f6f01f7940cdabf3afa4d351b7bb

作者:Piero Paialunga

目录
相关文章
|
3天前
|
缓存 开发者 Python
探索Python中的装饰器:从入门到实践
【9月更文挑战第36天】装饰器,在Python中是一种特殊的语法糖,它允许你在不修改原有函数代码的情况下,增加额外的功能。本文将通过浅显易懂的语言和实际代码示例,带你了解装饰器的基本原理,探索其背后的魔法,并展示如何在实际项目中运用这一强大工具。无论你是初学者还是有一定经验的开发者,这篇文章都将为你打开一扇通往更高效、更优雅代码的大门。
24 11
|
3天前
|
Python
Python 脚本高级编程:从基础到实践
本文介绍了Python脚本的高级概念与示例,涵盖函数的灵活应用、异常处理技巧、装饰器的使用方法、上下文管理器的实现以及并发与并行编程技术,展示了Python在自动化任务和数据操作中的强大功能。包括复杂函数参数处理、自定义装饰器、上下文管理器及多线程执行示例。
28 5
|
1天前
|
数据采集 数据可视化 数据挖掘
使用Python进行数据分析:从入门到实践
使用Python进行数据分析:从入门到实践
8 2
|
1天前
|
缓存 测试技术 开发者
探索Python中的装饰器:提升代码的灵活性和可维护性
在Python编程中,装饰器是一种强大且灵活的工具,它允许开发者在不修改现有代码的基础上,为函数或类添加额外的功能。本文将深入探讨装饰器的定义、使用场景以及如何创建自定义装饰器。通过实用的示例,我们将展示如何利用装饰器来增强代码的可重用性和可读性。
|
1天前
|
小程序 iOS开发 MacOS
将Python代码转化为可执行的程序
将Python代码转化为可执行的程序
10 1
|
1天前
|
存储 人工智能 开发工具
AI助理化繁为简,速取代码参数——使用python SDK 处理OSS存储的图片
只需要通过向AI助理提问的方式输入您的需求,即可瞬间获得核心流程代码及参数,缩短学习路径、提升开发效率。
52 0
AI助理化繁为简,速取代码参数——使用python SDK 处理OSS存储的图片
|
2天前
|
缓存 数据安全/隐私保护 开发者
探索Python中的装饰器:提升代码效率与可读性
在本文中,我们将深入探讨Python中的装饰器——一种能够修改或增强函数行为的强大工具。通过详细讲解装饰器的定义、使用方法以及实际案例分析,我们希望能够帮助读者更好地理解和应用这一技术,从而编写出更加高效和易读的代码。无论是初学者还是经验丰富的开发者,都能从本文中获得有价值的见解和技巧。
|
1天前
|
算法 开发者 计算机视觉
燃爆全场!Python并查集:数据结构界的网红,让你的代码炫酷无比!
在编程的世界里,总有一些数据结构以其独特的魅力和高效的性能脱颖而出,成为众多开发者追捧的“网红”。今天,我们要介绍的这位明星,就是Python中的并查集(Union-Find)——它不仅在解决特定问题上大放异彩,更以其优雅的设计和强大的功能,让你的代码炫酷无比,燃爆全场!
7 0
|
1天前
|
调度 Python
探索Python中的异步编程:从入门到实践
【8月更文挑战第70天】在Python的世界中,异步编程是一个能够显著提高程序性能和响应能力的技术。本文将通过一个简单的例子,介绍如何在Python中实现异步编程,以及如何利用这一技术优化你的代码。我们将从基础概念出发,逐步深入到实战应用,让你轻松掌握Python异步编程的精髓。
|
3天前
|
数据可视化 Python
Python 高级绘图:从基础到进阶的可视化实践
本文介绍了使用 Python 的强大绘图库 matplotlib 实现多种图表绘制的方法,包括简单的折线图、多条折线图、柱状图、饼图、散点图及 3D 图的绘制。通过具体代码示例展示了如何设置轴标签、标题、图例等元素,并指出了 matplotlib 支持更多高级绘图功能。来源:https://www.wodianping.com/app/2024-10/47112.html。
29 0