Numpy.fft.fft的功率谱

Question

无论我如何更改数据，通过下面的代码绘制的图形只是在零点周围波动。我的数据只有一列，记录了某种信号的每个时间点。 time_step 是我应该根据数据中两个相邻点的间隔定义的值吗？

data=np.loadtxt("timesequence",delimiter=",",usecols=(0,),unpack=True)

ps = np.abs(np.fft.fft(data))**2
time_step = 1

freqs = np.fft.fftfreq(data.size, time_step)
idx   = np.argsort(freqs)

pl.plot(freqs[idx], ps[idx])
pl.show()

- questionhang

---

我提供给 numpy.fft 的数据不是我接收信号的时间点？ - questionhang

这意味着你的 _f(t) = t_，它的傅里叶变换是 dirac delta 的 _一阶导数_。如果你在每个时间步骤接收一个信号，那么 data = [1, 1, 1, 1] 就是你的 信号 在 fft 中应该的值。 - askewchan

@questionhang 不，我认为那是你的问题。fft需要输入信号，你可以使用fftfreq将时间点转换为频率轴，以在功率谱图上显示。我已经为你提供了一个示例来演示这一点。 - Hooked

我创建了两列数据。一列像[1,1,1,1,...,1]，另一列像[1,2,3,4,5]。结果是一样的，在零附近有一个峰值。 - questionhang

如果你的信号在0附近不是大致对称的，那么具有高直流分量（fft的索引0）是正常的。尝试通过减去所有样本的平均值来预处理信号，以便你拥有大约相同数量的正能量和负能量... - twalberg

Answer 1

如其他人所提示，你的信号必须具有大的非零成分。在0（DC）处出现的峰值表示信号的平均值。这是由傅里叶变换本身导出的。余弦函数cos(0)*ps(0)表示信号的平均值的度量。傅里叶变换的其他组成部分是振幅不同的余弦波，它们显示了这些值的频率内容。

请注意，静止信号不会有很大的DC分量，因为它们已经是零均值信号。如果您不想有一个大的DC分量，则应计算信号的平均值并从中减去值。无论您的数据是0，...，999还是1，...，1000，甚至是1000，...，2000，您都将在0Hz处获得峰值。唯一的区别是峰值的大小，因为它测量的是平均值。

data1 = arange(1000)
data2 = arange(1000)+1000
dataTransformed3 = data - mean(data)
data4 = numpy.zeros(1000)
data4[::10] = 1 #simulate a photon counter where a 1 indicates a photon came in at time indexed by array. 
# we could assume that the sample rate was 10 Hz for example
ps1 = np.abs(np.fft.fft(data))**2
ps2 = np.abs(np.fft.fft(data))**2
ps3 = np.abs(np.fft.fft(dataTransformed))**2

figure()
plot(ps1) #shows the peak at 0 Hz
figure()
plot(ps2) #shows the peak at 0 Hz
figure()
plot(ps3) #shows the peak at 1 Hz this is because we removed the mean value but since
#the function is a step function the next largest component is the 1 Hz cosine wave.
#notice the order of magnitude difference in the two plots.

- Paul

---

我的数据是一系列光子的到达时间。如果单个光子均匀地到达，比如每秒钟一个计数，那么就存在一个周期，这个周期是一秒钟吗？如果光子不均匀地到达，我应该自己定义信号强度吗？ - questionhang

迷人的使用案例。 我认为无论哪种情况，你都没问题。 最后，您的数据数组可能不应该是arange（1000），而应该是类似于：data = zeros（1000）; data [:: 10] =1。这将表示每秒一次（假设采样率为10Hz-每10个值进入一个光子）会进入一个光子。 您可以直接测量与“理想”光子速率的偏差。 - Paul

我甚至不确定这是否是一个有意义的用例。听起来像是 OP 想要类似于 hist(diff(arrival_times)) 这样的东西，只是因为“我听说 FFT 与频率有关”。 - endolith

Answer 2

这是一个简单的示例，展示了输入和输出与峰值的关系，如您所期望的那样。

import numpy as np
from scipy.fftpack import rfft, irfft, fftfreq

time   = np.linspace(0,10,2000)
signal = np.cos(5*np.pi*time)

W = fftfreq(signal.size, d=time[1]-time[0])
f_signal = rfft(signal)

import pylab as plt
plt.subplot(121)
plt.plot(time,signal)
plt.subplot(122)
plt.plot(W,f_signal)
plt.xlim(0,10)
plt.show()

我使用rfft，因为很可能你的输入信号来自物理数据源，因此是实数。

- Hooked

---

是的，我的输入信号来自物理来源。绘制cos信号很容易。对于一系列的时间点，signal = np.cos(5np.pitime) 应该改为什么？ - questionhang

我每秒钟获得一个计数。它应该是周期性的。强度按计数/秒进行缩放，在这种完美的情况下，强度始终保持不变。 - questionhang

plt.plot(time,signal)和plt.plot(W,f_signal)都依赖于signal的定义。 - questionhang

如果signal是一个常数且W = fftfreq(signal.size, d=time[1]-time[0])，那么plt.plot(W,f_signal)就被确定了。这意味着时间点并不重要吗？ - questionhang

如果你使用subplot(122)，那么峰值应该在2pi/5pi=0.4左右。 - questionhang

Answer 3

如果您将数据全部变为正数：

ps = np.abs(np.fft.fft(data))**2
time_step = 1

如果您的实际数据振幅相对于大量的“DC”或0 Hz分量较小，那么很可能会创建一个大型的“DC”。因此，通过自动缩放功能，它将从图形中消失。

- jcoppens