Python - FFT 导致错误的物理含义答案

【问题标题】：Python - FFT leads to wrong physical meaningsPython - FFT 导致错误的物理含义
【发布时间】：2013-05-28 09:09:55
【问题描述】：

我是 Python 新手。我打算对一组离散点（时间、加速度）进行傅里叶变换，并将结果绘制出来。

我复制并粘贴示例 FFT 代码，并进行相应修改。

请看代码：

import numpy as np
import matplotlib.pyplot as plt

# Load the .txt file in
myData = np.loadtxt('twenty_z_up.txt')

# Extract the time and acceleration columns
time = copy(myData[:,0])

# Extract the acceleration columns
zAcc = copy(myData[:,3])

t = np.arange(10080)
sp = np.fft.fft(zAcc)
freq = np.fft.fftfreq(t.shape[-1])
plt.plot(freq, sp.real)

myData 是一个 10080 行 10 列的矩形矩阵。

因此，zAcc 是从矩阵中提取的第 3 行。

在 Spyder 绘制的图中，大部分谐波集中在 0 附近。它们都非常小。

但我的数据实际上是步行者携带手机的加速度（包括重力）。所以我预计最重要的谐波发生在 2Hz 附近。

为什么图表是无意义的？

提前致谢！

==============更新：我的图表======================

第一次域一：

x 轴以毫秒为单位。

y 轴以 m/s^2 为单位，由于地球重力，它的 DC 偏移量约为 10。

【问题讨论】：

np.arange(256) 创建一个整数数组（例如，从 0 到 255）。当您调用t.shape[-1] 时，您将获得数组t 的长度（即256）。 np.fft 是 python 快速傅立叶变换模块。 np.fft.fft() 是模块中的快速傅立叶变换函数。您需要在数据上调用该函数。
@nicholaschris 太棒了！谢谢！然后我将 256 修改为 10080 以适应我的情况。可以在更改后绘制该图。但是图表没有意义..你能帮我再检查一下吗？
0 附近的谐波是什么意思？如果您的意思是在 0Hz 处有能量，即 DC，这表示时域数据中存在 DC 偏移。上传图片会很有用（我不确定需要多少分才能将图片添加到问题中，如果您没有足够的链接，请留下指向 Imgur 或类似内容的链接，有人会为您编辑它）。
您并没有真正提供足够的信息来获得适当的帮助。毫无疑问，有人可以猜到这个问题，但这并不是微不足道的。时域是什么样的？如果您认为数据有明显的 2Hz 振荡，那在该图上应该很明显。
大部分谐波集中在 0 附近。它们都非常小。 - 你的轴刻度是多少？ 2Hz 可以接近 0Hz，尤其是在对数尺度上。

标签： python numpy python-3.x matplotlib

【解决方案1】：

您确实在（大约）2Hz 处获得了两个尖峰。您的采样周期约为 2.8 毫秒（我可以从您的第一个图中推断出最好的结果），给 +/-2Hz 的归一化频率为 +/-0.056，这大约是您的尖峰所在的位置。 fft.fftfreq 默认返回归一化频率（缩放采样周期）。您可以将d 参数设置为采样周期，您将得到一个包含实际频率的向量。

中间的巨大尖峰显然是 DC 偏移量（您可以通过减去平均值轻松消除）。

【讨论】：

哇。惊人的！是的。确实。但是由于我是 Python 新手，能否请您告诉我可以使用哪些 Python 语句来正确显示图表？通过“正确显示”，我的意思是 1. 删除琐碎的 DC； 2. 使 x 轴真正以 Hz 比例而不是 kHz 比例...
zAcc -= np.mean(zAcc) 去除平均值（使用就地操作），freq = np.fft.fftfreq(t.shape[-1], d=2.8e-3) 去除频率（以赫兹为单位）（d 是正确的采样周期）。此外，您使用副本，但不清楚它来自哪里。我建议你使用复制方法而不是复制功能：time = myData[:,0].copy()
太棒了！现在我明白了。非常感谢您的热心帮助！
另外，一切都很清楚，您的 FFT 图可能是 fft 的实部。我认为该情节放弃了虚构的部分，因为它无法合理地绘制它。

【解决方案2】：

正如其他人所说，我们需要查看数据，并将其发布到某个地方。只是为了检查，尝试首先在 fftfreq 中固定时间步长，然后绘制这个合成信号，然后绘制你的信号以查看它们的比较：

timestep=1./50.#Assume sampling at 50Hz. Change this accordingly.
N=10080#the number of samples
T=N*timestep
t = np.linspace(0,T,N)#needed only to generate xAcc_synthetic
freq=2.#peak a frequency at 2Hz
#generate synthetic signal at 2Hz and add some noise to it
xAcc_synthetic = sin((2*np.pi)*freq*t)+np.random.rand(N)*0.2
sp_synthetic = np.fft.fft(xAcc_synthetic)
freq = np.fft.fftfreq(t.size,d=timestep)
print max(abs(freq))==(1/timestep)/2.#simple check highest freq.
plt.plot(freq, abs(sp_synthetic))
xlabel('Hz')

现在，在等于 2 的 x 轴上，您实际上有一个 2Hz 的物理频率，您可能会发现您正在寻找的更明显的峰值。此外，您可能还想看看 yAcc 和 zAcc。

【讨论】：

我已经发布了比较图。请帮忙。 :)