从离散信号计算 FFT 特征

Question

我有一组以 10Hz 采样的均匀信号（意味着两个连续的数据点相隔 100 毫秒）。这实际上是 3d 陀螺仪的 3 个轴的大小，阵列包含 30 个数据点（在 3 秒内）。 我将这个系列的频率绘制如下

import numpy as np
import matplotlib.pyplot as pl

sample_rate = 10
x = np.array([318.45,302.78,316.47,334.14,333.41,326.15,320.07,318.68,314.12,308.64,300.15,304.33,318.42,322.72,329.56,339.18,338.03,343.27,351.44,353.23,352.35,352.88,353.43,352.14,351.28,352.82,353.36,353.35,353.19,353.82])
x = np.array(x) - np.mean(x)
p = np.abs(np.fft.rfft(x))
f = np.linspace(0, sample_rate/2, len(p))
pl.plot(f, p)
pl.show()

谁能告诉我我的设计是否正确？ 我打算计算以下特征（从上面的信号）

1. 直流分量
2. 光谱能量
3. 信息熵
4. 主要频率成分
5. 主要频率
6. FFT 分析前五个分量的大小

有人可以帮我填写上面的代码来计算这些特征吗？

----@RoadRunner66：请看下面我的问题，因为我无法给你回复很长的回复----

感谢您的回答和您的代码，

关于您的问题，数据来自测量欧拉角的陀螺仪。

那么 (sum x[i]**2 : 3357757.0) 是光谱能量吗？如果是，那么我是否需要通过将此数字除以 n 来对其进行标准化？ （或像你那样乘以 n），但是下面两篇论文的定义有所不同。

正如在第一篇论文（下面的第一个链接）中所说，“第二个频域特征集被选择为频谱能量，它被定义为平方 FFT 系数的总和”

在第二篇论文（第二个链接）中，他们以另一种方式表示“光谱能量：光谱系数的平方和除以窗口中的样本数”

那么主频率呢，它与主频率的含义（术语）相同吗？我猜主频率是指唯一具有最高频谱峰值的频率？

我像这样将频率和等效幅度打印成两行

我认为您像黄色波纹管一样打印了前 5 的大小。我不确定“前 5 个组件”的定义

如果我们像您指出的那样使用前连续的五个，那么包含那些（例如频率 0 或 0.666）并将它们输入我的预测模型（如下所述）是否有意义，因为它太低了与其他人相比。如果返回频谱清晰，主频为 1hz 和 3hz，那么 0.5hz 或 1.5hz 频率的幅度可能接近于零。

难道“FFT 分析的前五个分量的幅度”是我用蓝色突出显示的“前五个主要分量的幅度”吗？这个术语是指 5 个值还是仅指 1 个值（5 个平方和的平方根）？ 如果它指的是 5 个值（很可能是），那么在比较两个信号之间的差异时，我认为幅度上排名前 5 位的主要分量将是一个更好的选择？

顺便说一句，第二篇论文也写了“前5-FFT系数：前5个快速傅里叶变换系数是因为它们捕获了主要的频率分量，并且使用附加系数并没有提高准确度”

坦率地说，我正在处理奶牛活动的问题，我的策略是将传感器数据分割成时间窗口（3、5、7..秒）并从每个窗口中提取特征，然后将它们提供给机器学习模型。

（我的数据包括一个 3d 陀螺仪和一个附在牛脖子上的 3d 加速度计，传感器数据采样为 10Hz）

我想结合两种特征，一种是时域特征，一种是频域特征。

我阅读了论文并发现了上述一组频域特征，其中包括术语“FFT 分析的前五个分量的幅度”（来自本文https://ieeexplore.ieee.org/abstract/document/4663615） 从这个https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4634510/

他们将第二个称为“第一个 5-FFT 系数：前 5 个快速傅里叶变换系数被采用，因为它们捕获了主要的频率分量，并且使用附加系数并没有提高准确性。”

非常感谢您的阅读和回答！

Answer 1

您已经完成了大部分工作。这似乎是针对班级分配的，因此这可能解释了数据的稀疏性，而这些数据实际上并没有太多可取之处。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

sample_rate = 10 #Hz

x = np.array([318.45,302.78,316.47,334.14,333.41,326.15,320.07,318.68,314.12,308.64,
              300.15,304.33,318.42,322.72,329.56,339.18,338.03,343.27,351.44,353.23,
              352.35,352.88,353.43,352.14,351.28,352.82,353.36,353.35,353.19,353.82])
n=len(x)
print(f'number of points n : {n}')
t= np.linspace(0,3, n)     # time according to your description
mn=np.mean(x)
print(f' mean = {mn:.3f} unit')
print(f' sum x[i]**2  : {np.sum(x**2) :.1f} unit^2 ')

pr = np.abs(np.fft.rfft(x))/n
print()
print(f' DC peak :  {pr[0]};  this makes choice of normalization 1/n meaningful')

print(f' n *sum X[k]**2   : {np.sum(np.abs(pr)**2)*n :.1f} unit^2 ') 
f = np.linspace(0, sample_rate/2, len(pr))   # 10 Sa/sec, so 5 Hz is Nyquist limit

plt.figure(figsize=(20,5))

plt.subplot(141)
plt.plot(t,x,'.-')
plt.title('original data') 

plt.subplot(142)
plt.plot(f, pr,'.-')
plt.title('spectrum')


plt.subplot(143)
xf = np.array(x) - mn   #  remove the DC
plt.plot(t, xf,'.-')
plt.title('original data with DC removed')

plt.subplot(144)
pr = np.abs(np.fft.rfft(xf))/n
plt.plot(f, pr, '.-');
plt.title('spectrum (DC removed) ')

直流峰值与原始信号在 334 个单位处的平均值相同。

频谱能量（参见Parseval's Theorem）在时域和频域中是相同的。

我不知道您应该使用什么信息熵定义（参见例如Approx. entropy）。

频谱中的主要频率（去除 DC 后）为 1/3 Hz、1 Hz 和 2 Hz，其中 1/3 Hz 最大。我打印了前五个的大小。

对我来说一个重要的问题是数据的物理意义。它们是角度吗？如果有，在哪个单位？

请注意，可以“看到”原始数据中 1/3 Hz（一波超过 3 秒）、1 Hz（一波超过 1 秒）和 2 Hz（一波超过 1/2 秒）的频率分量)，分别。