看来 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)
Now 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） 数据。给出的频率为lombscargle不等于1/trange, where trange是数据的完整时间跨度，1/dt.min()。频率数为16000（num_ls_freqs在脚本中）。一些试验表明，这大致是解析频谱所需的最小数量。小于这个值，山峰就会开始移动。除此之外，频谱几乎没有变化。计算表明在0.0242 Hz处有一个峰值。其他峰值是该频率的谐波。
现在我们已经估计了基频，我们可以用它来指导直方图中 bin 大小的选择timepoints用于 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是从文件中读取的原始时间集；它包括许多重复的时间值，如上所述。）
然后我们应用FFThist。 （我们实际上会使用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()
