有人能指出为什么我得到的结果截然不同吗？
有很多不应该出现的峰。事实上，应该只有一个峰。
我是一个 python 新手，欢迎对下面我的代码提出所有评论。
测试数据在这里。 enter link description here
您可以直接 wget https://clbin.com/YJkwr 。
它的第一列是一系列光子的到达时间。光源随机发射光子。
总时间为55736s，有67398个光子。
我尝试检测光强度的某种周期性。
我们可以对时间进行分箱，光强度与每个时间分箱中的光子数成正比。
我尝试了 scipy.signal 的 numpy.fft 和 lomb-scargle 来制作其功率谱，但得到了截然不同的结果。
快速傅立叶变换
import pylab as pl
import numpy as np
timepoints=np.loadtxt('timesequence',usecols=(0,),unpack=True,delimiter=",")
binshu=50000
interd=54425./binshu
t=np.histogram(timepoints,bins=binshu)[0]
sp = np.fft.fft(t)
freq = np.fft.fftfreq(len(t),d=interd)
freqnum = np.fft.fftfreq(len(t),d=interd).argsort()
pl.xlabel("frequency(Hz)")
pl.plot(freq[freqnum],np.abs(sp)[freqnum])
lomb-scargle
timepoints=np.loadtxt('timesequence',usecols=(0,),unpack=True,delimiter=",")
binshu=50000
intensity=np.histogram(timepoints,bins=binshu)[0].astype('float64')
middletime=np.histogram(timepoints,bins=binshu)[1][1:binshu+1]-np.histogram(timepoints,bins=binshu)[1][3]*0.5
freq1=1./(timepoints.max()-timepoints.min())
freq2=freq1*len(timepoints)
freqs=np.linspace(freq1,freq2,1000)
pm1=spectral.lombscargle(middletime,intensity,freqs)
pl.xlabel("frequency(Hz)")
pl.plot(freqs,pm1)
pl.xlim(0.5*freq1,2*freq2)
pl.ylim(0,250)
pl.show()
************************************
谢谢你，沃伦，我很感激。谢谢你的详细回复。
你说得对，有一个积分时间，这里大约是1.7s。
单次曝光次数较多，每次曝光耗时1.7s
在一次曝光中，我们无法准确判断其到达时间。
如果时间序列是这样的:
0 1.7 3.4 8.5 8.5
最后两个光子的积分时间是1.7s，而不是(8.5-3.4)s。所以我会修改你的部分代码。
但是我的问题仍然存在。您可以调整几个参数，以在某种程度上获得 lomb-scargle 峰值中的 0.024Hz 峰值。您可以使用它来指导 fft 中的参数。
如果您不知道数字0.024，也许您可以使用不同的参数来获得不同的最高峰？
如何保证每次都能得到正确的num_ls_freqs？您可以看到如果我们选择不同的 num_ls_freqs，最高峰值偏移。
如果我有很多计时系列，每次我都必须指定不同的参数？以及如何获得它们？
最佳答案
看来 50000 个 bin 还不够。如果你改变binshu ，您会看到尖峰的位置发生了变化，而且变化不只是一点点。
在深入研究 FFT 或 Lomb-Scargle 之前，我发现首先研究原始数据很有用。以下是文件的开头和结尾:
0,1,3.77903
0,1,4.96859
0,1,1.69098
1.74101,1,4.87652
1.74101,1,5.15564
1.74101,1,2.73634
3.48202,1,3.18583
3.48202,1,4.0806
5.2229,1,1.86738
6.96394,1,7.27398
6.96394,1,3.59345
8.70496,1,4.13443
8.70496,1,2.97584
...
55731.7,1,5.74469
55731.7,1,8.24042
55733.5,1,4.43419
55733.5,1,5.02874
55735.2,1,3.94129
55735.2,1,3.54618
55736.9,1,3.99042
55736.9,1,5.6754
55736.9,1,7.22691
存在相同到达时间的集群。 (这是光子探测器的原始输出，还是该文件已以某种方式进行了预处理？)
因此第一步是将这些簇一次性合并为一个数字，因此数据看起来像(时间点，计数)，例如:
0, 3
1.74101, 3
3.48202, 2
5.2229, 1
etc
以下将完成计数:
times, inv = np.unique(timepoints, return_inverse=True)
counts = np.bincount(inv)
现在times是一个递增的唯一时间值的数组，并且 counts是当时检测到的光子数。 (请注意，我猜测这是对数据的正确解释!)
让我们看看采样是否接近均匀。到达间隔时间数组为
dt = np.diff(times)
如果采样完全均匀，则 dt 中的所有值会是一样的。我们可以在timepoints上使用上面使用的相同模式。查找 dt 中的唯一值(及其出现频率) :
dts, inv = np.unique(dt, return_inverse=True)
dt_counts = np.bincount(inv)
例如，这就是 dts如果我们打印它看起来像:
[ 1.7 1.7 1.7 1.7 1.74 1.74 1.74 1.74
1.74 1.74 1.74088 1.741 1.741 1.741 1.741
1.741 1.741 1.741 1.74101 1.74102 1.74104 1.7411
1.7411 1.7411 1.7411 1.7411 1.742 1.742 1.742
1.746 1.75 1.8 1.8 1.8 1.8 3.4 3.4
3.4 3.4 3.48 3.48 3.48 3.48 3.48 3.482
3.482 3.482 3.482 3.482 3.4821 3.483 3.483 3.49
3.49 3.49 3.49 3.49 3.5 3.5 5.2 5.2
5.2 5.2 5.22 5.22 5.22 5.22 5.22 5.223
5.223 5.223 5.223 5.223 5.223 5.23 5.23 5.23
5.3 5.3 5.3 5.3 6.9 6.9 6.9 6.9
6.96 6.96 6.964 6.964 6.9641 7. 8.7 8.7
8.7 8.7 8.7 8.71 10.4 10.5 12.2 ]
(表面上的重复值实际上是不同的。未显示数字的完整精度。)如果采样完全均匀，则该列表中将只有一个数字。相反，我们看到的是数字簇，没有明显的主导值。 (1.74 的倍数占优势——这是探测器的特性吗？)
根据这一观察，我们将从 Lomb-Scargle 开始。下面的脚本包含用于计算和绘制 ( times , counts ) 数据的 Lomb-Scargle 周期图的代码。 lombscargle 的频率不同于 1/trange ，其中trange是数据的完整时间跨度，到 1/dt.min() 。频率数为 16000(脚本中的 num_ls_freqs)。一些试验表明，这大致是解析频谱所需的最小数量。小于这个值，山峰就会开始移动。除此之外，频谱几乎没有变化。计算表明在0.0242 Hz处有一个峰值。其他峰值是该频率的谐波。
现在我们已经估计了基频，我们可以用它来指导 timepoints 直方图中 bin 大小的选择用于 FFT 计算。我们将使用导致基频过采样 8 倍的 bin 大小。(在下面的脚本中，这是 m 。)也就是说，我们这样做
m = 8
nbins = int(m * ls_peak_freq * trange + 0.5)
hist, bin_edges = np.histogram(timepoints, bins=nbins, density=True)
( timepoints 是从文件中读取的原始时间集；它包括许多重复的时间值，如上所述。)
然后我们将 FFT 应用于 hist 。 (我们实际上将使用 numpy.fft.rfft 。)以下是使用 FFT 计算的频谱图:
正如预期的那样，在 0.0242 Hz 处有一个峰值。
这是脚本:
import numpy as np
from scipy.signal import lombscargle
import matplotlib.pyplot as plt
timepoints = np.loadtxt('timesequence', usecols=(0,), delimiter=",")
# Coalesce the repeated times into the `times` and `counts` arrays.
times, inv = np.unique(timepoints, return_inverse=True)
counts = np.bincount(inv)
# Check the sample spacing--is the sampling uniform?
dt = np.diff(times)
dts, inv = np.unique(dt, return_inverse=True)
dt_counts = np.bincount(inv)
print dts
# Inspection of `dts` shows that the smallest dt is about 1.7, and there
# are many near multiples of 1.74, but the sampling is not uniform,
# so we'll analyze the spectrum using lombscargle.
# First remove the mean from the data. This is not essential; it just
# removes the large value at the 0 frequency that we don't care about.
counts0 = counts - counts.mean()
# Minimum interarrival time.
dt_min = dt.min()
# Total time span.
trange = times[-1] - times[0]
# --- Lomb-Scargle calculation ---
num_ls_freqs = 16000
ls_min_freq = 1.0 / trange
ls_max_freq = 1.0 / dt_min
freqs = np.linspace(ls_min_freq, ls_max_freq, num_ls_freqs)
ls_pgram = lombscargle(times, counts0, 2*np.pi*freqs)
ls_peak_k = ls_pgram.argmax()
ls_peak_freq = freqs[ls_peak_k]
print "ls_peak_freq =", ls_peak_freq
# --- FFT calculation of the binned data ---
# Assume the Lomb-Scargle calculation gave a good estimate
# of the fundamental frequency. Use a bin size for the histogram
# of timepoints that oversamples that period by m.
m = 8
nbins = int(m * ls_peak_freq * trange + 0.5)
hist, bin_edges = np.histogram(timepoints, bins=nbins, density=True)
delta = bin_edges[1] - bin_edges[0]
fft_coeffs = np.fft.rfft(hist - hist.mean())
fft_freqs = np.fft.fftfreq(hist.size, d=delta)[:fft_coeffs.size]
# Hack to handle the case when hist.size is even. `fftfreq` puts
# -nyquist where we want +nyquist.
fft_freqs[-1] = abs(fft_freqs[-1])
fft_peak_k = np.abs(fft_coeffs).argmax()
fft_peak_freq = fft_freqs[fft_peak_k]
print "fft_peak_freq =", fft_peak_freq
# --- Lomb-Scargle plot ---
plt.figure(1)
plt.clf()
plt.plot(freqs, ls_pgram)
plt.title('Spectrum computed by Lomb-Scargle')
plt.annotate("%6.4f Hz" % ls_peak_freq,
xy=(ls_peak_freq, ls_pgram[ls_peak_k]),
xytext=(10, -10), textcoords='offset points')
plt.xlabel('Frequency (Hz)')
plt.grid(True)
# --- FFT plot ---
plt.figure(2)
plt.clf()
plt.plot(fft_freqs, np.abs(fft_coeffs)**2)
plt.annotate("%6.4f Hz" % fft_peak_freq,
xy=(fft_peak_freq, np.abs(fft_coeffs[fft_peak_k])**2),
xytext=(10, -10), textcoords='offset points')
plt.title("Spectrum computed by FFT")
plt.grid(True)
plt.show()
关于numpy - 我的功率谱可信吗？ lomb-scargle 和 fft 之间的比较(scipy.signal 和 numpy.fft)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/20658318/