如何计算给定波的频率和时间

Question

我有速度与时间的数据。时间步长不一致，但是力度数据是波动的。如何使用Python的FFT计算速度的主频率？我在网上看到的大多数示例都是为了统一时间。

我的数据就像

7.56683038E+02  2.12072850E-01 
7.56703750E+02  2.13280844E-01
7.56724461E+02  2.14506402E-01
7.56745172E+02  2.15748934E-01
7.56765884E+02  2.17007907E-01
7.56786595E+02  2.18282753E-01

像这样的10000行。

看到一些在线响应，我编写了如下代码，但它不起作用：

#!/usr/bin/env python

import numpy as np
import scipy as sy
import scipy.fftpack as syfp
import pylab as pyl

# Calculate the number of data points
length = 0
for line in open("data.dat"):
    length = length + 1

# Read in data from file here
t = np.zeros(shape=(length,1))
u = np.zeros(shape=(length,1))


length = 0
for line in open("data.dat"):
    columns = line.split(' ')
    t[length] = float(columns[0])
    u[length] = float(columns[1])
    length = length + 1

# Do FFT analysis of array
FFT = sy.fft(u)

# Getting the related frequencies
freqs = syfp.fftfreq(len(u))

# Create subplot windows and show plot
pyl.subplot(211)
pyl.plot(t, u)
pyl.xlabel('Time')
pyl.ylabel('Amplitude')
pyl.subplot(212)
pyl.plot(freqs, sy.log10(FFT), 'x')
pyl.show()

----------------------编辑------------------------

通过此代码，我得到如下图所示的输出。我不确定这个数字显示了什么。我期望在FFT图中看到一个峰值

----------------------编辑------------------------

我在以下注释中建议的带有sin函数的模拟数据结果如下：

Answer 1

从我所看到的来看，您的代码基本上不错，但是缺少一些细节。我认为您的问题主要与解释有关。因此，现在最好看一下模拟数据，这是我在注释中建议的模拟数据示例（并且我已经添加了有关重要行和##更改的注释）：

import numpy as np
import scipy as sy
import scipy.fftpack as syfp
import pylab as pyl

dt = 0.02071 
t = dt*np.arange(100000)             ## time at the same spacing of your data
u = np.sin(2*np.pi*t*.01)            ## Note freq=0.01 Hz

# Do FFT analysis of array
FFT = sy.fft(u)

# Getting the related frequencies
freqs = syfp.fftfreq(len(u), dt)     ## added dt, so x-axis is in meaningful units

# Create subplot windows and show plot
pyl.subplot(211)
pyl.plot(t, u)
pyl.xlabel('Time')
pyl.ylabel('Amplitude')
pyl.subplot(212)
pyl.plot(freqs, sy.log10(abs(FFT)), '.')  ## it's important to have the abs here
pyl.xlim(-.05, .05)                       ## zoom into see what's going on at the peak
pyl.show()

如您所见，在输入频率（.01 Hz）的+和-处有两个峰值，符合预期。

编辑：
困惑为什么这种方法不适用于OP的数据，我也对此进行了研究。问题在于采样时间间隔不均匀。这是时间的直方图（下面的代码）。

因此，样本之间的时间大致在短时间和长时间之间平均分配。我快速浏览了一下这里的模式，没有什么是显而易见的。

要进行FFT，需要均匀间隔的时间样本，因此我进行插值得到以下内容：

这是合理的（直流偏移，初级峰值和小谐波）。这是代码：

data = np.loadtxt("data.dat", usecols=(0,1))
t = data[:,0]
u = data[:,1]

tdiff = t[1:]-t[:-1]
plt.hist(tdiff)

new_times = np.linspace(min(t), max(t), len(t))
new_data = np.interp(new_times, t, u)

t, u = new_times, new_data

FFT = sy.fft(u)
freqs = syfp.fftfreq(len(u), dt)

# then plot as above