从离散信号计算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()

有人可以告诉我我的计划正确吗？我正计划计算以下特征（从上方信号）
直流分量
光谱能量
信息熵
主频成分
主频率
FFT分析的前五个分量的幅度

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

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) ')

DC峰值与原始信号在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秒内一波）“看到”原始数据中的频率分量分别。