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

Posted 2023-03-12

【中文标题】如何计算给定波的频率和时间

【英文标题】：How to calculate frequency of a give wave and time

【发布时间】：2014-11-24 04:31:08

【问题描述】：

我有速度与时间的数据。时间步长不均匀，但速度数据是一个波。如何使用 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 图中看到一个峰值

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

我的模拟数据与下面 cmets 中建议的 sin 函数的结果在这里：

【问题讨论】：

“它不工作”是什么意思？另外，“主频”是什么意思？

在FFT图中，如果你绘制sy.log10(np.abs(FFT))，你会看到什么（即，注意那里abs的使用）。

这不会改变输出中的任何内容。顺便说一句，我检查了一下，时间步长看起来几乎一致，dt= 0.02071，我的文件中有 100000 个数据点。

很难用我们没有的数据来诊断结果。尝试一些假数据怎么样，你知道结果应该是什么样子：t = 0.02071*np.arange(100000) 和u = np.sin(2*np.pi*t*.01)？

感谢您的想法。我会试试的。同时，我在drive.google.com/file/d/0B9w1Yb4SFBxrWEFGTERVX1dZSGM/… 上提供了数据。第一列是时间，第二列是速度。不要打扰其他列。

【参考方案1】：

据我所知，您的代码基本没问题，但缺少一些细节。我认为您的问题主要与解释有关。正因为如此，mock 数据是现在最好看的，这里有一个我在 cmets 中建议的 mock 数据的例子（我已经添加了关于重要行的 cmets，## 进行更改）：

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

【讨论】：

感谢您的回复。看起来代码确实在做正确的事情。我做了你推荐的改变。但是，在我输入数据之后，我得到的 FFT 的输出就像一个波（如我的问题的第一个图所示）。我的原始数据看起来像一个平滑的波，所以我不知道如何解释我的输出。

@jhaprade：我最终对此感到困惑，所以我做到了，并在编辑中包含了您的数据问题的答案。