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

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图中看到一个峰值

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

我与在下面的意见建议的罪恶功能模拟数据结果在这里：

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() 

正如你可以看到，有两个峰，在+和 - 输入频率（0.01赫兹），如所预期。

编辑：
疑惑，为什么这种方法没有为OP的数据工作，我看了一下这一点。问题是采样时间不均匀。这是时间的直方图（下面的代码）。

So样本之间的时间被大致均匀地在短时间和长时间之间拆分。我在这里快速查找了一个模式，没有什么是显而易见的。

要执行的FFT，需要均匀间隔的时间样本，所以内插，以得到以下：

这是合理的（DC偏移，主峰值和小谐波）。这里的代码：

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