在Python中确定数组频率

Question

我有一个示例文件，里面全是浮点数，内容如下：

    -0.02  3.04  3.04  3.02  3.02  3.06  3.04  3.02  3.04  3.02  3.04  3.02
     3.04  3.02  3.04  3.04  3.04  3.02  3.04  3.02  3.04  3.02  3.04  3.02
     3.06  3.02  3.04  3.02  3.04  3.02  3.02  3.06  3.04  3.02  3.04  3.02
     3.04  3.02  3.04  3.04  3.04  3.02  3.04  3.02  3.02  3.06  3.04  3.02
     3.06  3.02  3.04 -0.02 -0.02 -0.02 -0.02 -0.02 -0.02 -0.04 -0.02 -0.04

这些数字被放在一个文本文件里。我想读取这个文本文件，并确定这个信号的频率。这些数据是从数字示波器上获取的。我可以在示波器的显示屏上看到频率，但我也想通过在Python中处理这些数据来验证一下。我是在电脑上用Python从设备获取数据的。

虽然我在Python中可以做一些基础的操作，但我对文本处理完全是个新手。我想我需要先把文件中的数据加载到一个数组里，然后再进行FFT（快速傅里叶变换）或者其他简单的算法，这样就能得到一个以赫兹为单位的整数频率。

理论上，我知道如何进行傅里叶分析，并且可以在纸上为任何给定的信号计算。但我不知道在Python中如何开始处理这些数据集。我已经尝试过查看scipy和numpy的文档，但对我来说效果不是很好。

我希望能得到一些有经验的用户的指导。

Answer 1

从你的问题中不太清楚文件里的值具体代表什么。不过假设这些值是连续的电压采样数据，你可以用下面的代码把文件加载到一个Numpy数组中：

import numpy as np
data = np.array([float(f) for f in file(filename).read().split()])

然后你可以用下面的代码计算傅里叶变换：

import numpy.fft as fft
spectrum = fft.fft(data)

接着，你可以用下面的代码绘制FFT的幅度：

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

你看到的结果应该和示波器上显示的内容一致。

如果你想找出频谱中主要的频率，你需要在某个阈值处切割数组，比如可以这样做：

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

freq（以及peaks）的内容是以采样率为单位的频率。例如，如果你的示波器每微秒采样一次波形，那么freq里的值就是以兆赫为单位的。所以如果你输入一个理想的1 kHz信号，比如可以用下面的代码：

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

你会在0.001 MHz的位置看到一个峰值，因此你会发现peaks = array([-0.001, 0.001])。