Python Numpy FFT - 或 - FFT 来查找波的周期而不是其频率？答案

【问题标题】：Python Numpy FFT -or- RFFT to find period of a wave instead of its frequiency?Python Numpy FFT - 或 - FFT 来查找波的周期而不是其频率？
【发布时间】：2014-12-28 03:37:13
【问题描述】：

我是信号分析的新手，我想参加一个项目，尝试通过分析我们实验室的气温稳定性来学习 Python 的 FFT 模块。

我编写了这个 Python 脚本，其中包含来自我们传感器的一些真实数据。我将在这里解释一些初始变量：

“数据” 是从数据库中获取的数据。通常可以假设它们以 120 秒为间隔，但不能保证。所以为了帮助计算我添加的快速平均采样率：

"temporal_window" 这是从第一次测量到最后一次测量的时间（以秒为单位）。那么在哪里：

T = temporal_window/N #should equal roughly 120 seconds

“调试” 在正常操作中，数据通过从数据库构建的数组（也称为“数据”）馈送到 FFT，但当我试图了解 FFT 的工作原理时，我决定制作一个“diagnostics_array”，它只是一个数组，其数据点数与数据库中的数组相同，但具有给定波长以秒为单位的正弦波。

import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt

data = np.array([17.38 , 17.66 , 18.26 , 18.62 , 18.98 , 19.42 , 19.7 , 19.38 , 18.46 , 17.82 , 17.5 , 17.3 , 17.9 , 18.3 , 18.66 , 19.06 , 19.5 , 19.78 , 19.94 , 19.06 , 18.06 , 17.54 , 17.26 , 18.02 , 18.42 , 18.78 , 19.18 , 19.54 , 19.82 , 19.42 , 18.54 , 17.74 , 17.34 , 17.18 , 17.86 , 18.38 , 18.7 , 19.02 , 19.42 , 19.7 , 19.42 , 18.38 , 17.74 , 17.34 , 17.66 , 18.22 , 18.46 , 18.82 , 19.26 , 19.62 , 19.78 , 18.78 , 17.98 , 17.46 , 17.3 , 17.98 , 18.38 , 18.74 , 19.06 , 19.42 , 19.74 , 19.98 , 19.54 , 18.46 , 17.82 , 17.26 , 17.7 , 18.3 , 18.62 , 18.98 , 19.42 , 19.74 , 19.9 , 19.1 , 18.14 , 17.74 , 17.98 , 18.38 , 18.74 , 19.1 , 19.54 , 19.82 , 19.38 , 18.54 , 17.9 , 17.58 , 18.14 , 18.58 , 18.9 , 19.3 , 19.62 , 19.9 , 19.54 , 18.54 , 17.82 , 17.38 , 17.74 , 18.3 , 18.7 , 19.1 , 19.42 , 19.66 , 18.78 , 17.94 , 17.42 , 17.22 , 17.94 , 18.38 , 18.82 , 19.18 , 19.58 , 19.82 , 19.94 , 19.02 , 18.22 , 17.66 , 17.46 , 18.1 , 18.46 , 18.86 , 19.18 , 19.58 , 19.9 , 19.46 , 18.5 , 17.82 , 17.38 , 17.66 , 18.26 , 18.66 , 19.02 , 19.46 , 19.78 , 19.94 , 19.06 , 19.18 , 19.58 , 19.94 , 20.22 , 20.38 , 20.54 , 20.58 , 20.06 , 18.94 , 18.14 , 17.74 , 17.34 , 17.7 , 18.3 , 18.7 , 19.02 , 19.42 , 19.74 , 19.9 , 19.02 , 18.22 , 17.66 , 17.3 , 17.7 , 18.3 , 18.7 , 18.98 , 19.38 , 19.74 , 19.42 , 18.5 , 17.74 , 17.26 , 17.66 , 18.3 , 18.62 , 19.02 , 19.42 , 19.74 , 19.94 , 18.98 , 18.22 , 17.78 , 17.58 , 18.14 , 18.5 , 18.86 , 19.18 , 19.58 , 19.78 , 18.86 , 18.02 , 17.58 , 17.34 , 18.02 , 18.38 , 18.78 , 19.14 , 19.58 , 19.82 , 19.5 , 18.5 , 17.86 , 17.46 , 17.74 , 18.3 , 18.62 , 19.06 , 19.42 , 19.74 , 18.86 , 17.98 , 17.54 , 17.18 , 17.98 , 18.38 , 18.74 , 19.1 , 19.54 , 19.86 , 19.46 , 18.46 , 17.9 , 17.3 , 17.66 , 18.22 , 18.66 , 18.94 , 19.42 , 19.78 , 19.42 , 18.46 , 17.82 , 18.02 , 18.5 , 18.86 , 19.26 , 19.62 , 19.34 , 18.42 , 17.86 , 18.02 , 18.46 , 18.78 , 19.26 , 19.58 , 19.34 , 18.3 , 17.7 , 17.42 , 18.1 , 18.5 , 18.78 , 19.22 , 19.62 , 19.74 , 18.78 , 17.98 , 17.42 , 17.14 , 17.42 , 18.02 , 18.42 , 18.74 , 19.14 , 19.5 , 19])
temporal_window = 42014.0 #seconds

N = len(data) #datapoints
T = temporal_window/N #should equal roughly 120 seconds

###Diagnostic Override###
debug = True #DEBUG SWITCH
if debug:
    wave_lenght = 60*60*1 #in seconds (eg. 60*60*2 = 2 hours)
    print "Created a sine wave with %s second period" % wave_lenght
    diagnostic_array = np.arange(0,1,1./N)
    diagnostic_array = np.cos(2*np.pi*temporal_window/wave_lenght*diagnostic_array)
    data = diagnostic_array
#########################

a=np.abs(fft.rfft(data))
a[0]=0 #Not sure if this is a good idea but seems to help with choppy data..
xt = np.linspace(0.0, temporal_window, a.size)

print "Peak found at %s second period" % int(xt[np.argmax(a)])

plt.subplot(211)
plt.plot(xt,a)
plt.subplot(212)
plt.plot(np.linspace(0,temporal_window,data.size),data)
plt.show()

所以当从上面运行代码时，我得到以下打印语句：

>>> #1 hour period
Created a sine wave with 3600 second period
Peak found at 3848 second period

>>> #2 hour period
Created a sine wave with 7200 second period
Peak found at 1924 second period

因此，随着波长变长（完全符合预期），FFT 峰值的结果似乎会变小。但我不确定如何改变它，以便在这个例子中峰值与波长相匹配，以秒为单位。 FFT可以吗？我正在阅读有关 IFFT 以转换回时域的信息，但对这个主题没有很好的理解，我有点不知所措..

任何关于如何实现这一点的想法或想法将不胜感激！如果我没有清楚地解释我的意图，请告诉我，我很乐意补充细节。非常感谢！！

【问题讨论】：

period = sample_rate / frequency，反之亦然frequency = sample_rate / period

标签： python numpy fft

【解决方案1】：

感谢 hobbs 的推动，我重新评估了我实际看到的内容。

经过进一步研究，我发现 rfftfreq 函数比 linspace 更方便。

因此，这是似乎按预期工作的更新代码。作为注释，我在执行 np.divide(60,freqs) 时得到“RuntimeWarning：除以零”。然而，这似乎不会影响结果。

我确实注意到，对于脚本的当前诊断部分，它允许 FFT 中的泄漏，因为它不关心将整个波拟合到数据集（例如，可能是 1.3 个波长或其他）。

所以要真正看到这一点（峰值 FFT 与输入波形周期匹配），您所要做的就是改变这一行：

-来自-

wave_lenght = 60*60*1 #in seconds (eg. 60*60*2 = 2 hours)

-to-

whole_waves = 2
wave_lenght = temporal_window/whole_waves #fits n number of whole waves within the dataset

这使得波成为总时间的函数，而不是设定波长，因此它很好地适合数据集。

这是完整的更新脚本。如果有人发现错误，请发表评论（我仍在学习这些东西并喜欢社区的反馈）！

import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt

data = np.array([17.38 , 17.66 , 18.26 , 18.62 , 18.98 , 19.42 , 19.7 , 19.38 , 18.46 , 17.82 , 17.5 , 17.3 , 17.9 , 18.3 , 18.66 , 19.06 , 19.5 , 19.78 , 19.94 , 19.06 , 18.06 , 17.54 , 17.26 , 18.02 , 18.42 , 18.78 , 19.18 , 19.54 , 19.82 , 19.42 , 18.54 , 17.74 , 17.34 , 17.18 , 17.86 , 18.38 , 18.7 , 19.02 , 19.42 , 19.7 , 19.42 , 18.38 , 17.74 , 17.34 , 17.66 , 18.22 , 18.46 , 18.82 , 19.26 , 19.62 , 19.78 , 18.78 , 17.98 , 17.46 , 17.3 , 17.98 , 18.38 , 18.74 , 19.06 , 19.42 , 19.74 , 19.98 , 19.54 , 18.46 , 17.82 , 17.26 , 17.7 , 18.3 , 18.62 , 18.98 , 19.42 , 19.74 , 19.9 , 19.1 , 18.14 , 17.74 , 17.98 , 18.38 , 18.74 , 19.1 , 19.54 , 19.82 , 19.38 , 18.54 , 17.9 , 17.58 , 18.14 , 18.58 , 18.9 , 19.3 , 19.62 , 19.9 , 19.54 , 18.54 , 17.82 , 17.38 , 17.74 , 18.3 , 18.7 , 19.1 , 19.42 , 19.66 , 18.78 , 17.94 , 17.42 , 17.22 , 17.94 , 18.38 , 18.82 , 19.18 , 19.58 , 19.82 , 19.94 , 19.02 , 18.22 , 17.66 , 17.46 , 18.1 , 18.46 , 18.86 , 19.18 , 19.58 , 19.9 , 19.46 , 18.5 , 17.82 , 17.38 , 17.66 , 18.26 , 18.66 , 19.02 , 19.46 , 19.78 , 19.94 , 19.06 , 19.18 , 19.58 , 19.94 , 20.22 , 20.38 , 20.54 , 20.58 , 20.06 , 18.94 , 18.14 , 17.74 , 17.34 , 17.7 , 18.3 , 18.7 , 19.02 , 19.42 , 19.74 , 19.9 , 19.02 , 18.22 , 17.66 , 17.3 , 17.7 , 18.3 , 18.7 , 18.98 , 19.38 , 19.74 , 19.42 , 18.5 , 17.74 , 17.26 , 17.66 , 18.3 , 18.62 , 19.02 , 19.42 , 19.74 , 19.94 , 18.98 , 18.22 , 17.78 , 17.58 , 18.14 , 18.5 , 18.86 , 19.18 , 19.58 , 19.78 , 18.86 , 18.02 , 17.58 , 17.34 , 18.02 , 18.38 , 18.78 , 19.14 , 19.58 , 19.82 , 19.5 , 18.5 , 17.86 , 17.46 , 17.74 , 18.3 , 18.62 , 19.06 , 19.42 , 19.74 , 18.86 , 17.98 , 17.54 , 17.18 , 17.98 , 18.38 , 18.74 , 19.1 , 19.54 , 19.86 , 19.46 , 18.46 , 17.9 , 17.3 , 17.66 , 18.22 , 18.66 , 18.94 , 19.42 , 19.78 , 19.42 , 18.46 , 17.82 , 18.02 , 18.5 , 18.86 , 19.26 , 19.62 , 19.34 , 18.42 , 17.86 , 18.02 , 18.46 , 18.78 , 19.26 , 19.58 , 19.34 , 18.3 , 17.7 , 17.42 , 18.1 , 18.5 , 18.78 , 19.22 , 19.62 , 19.74 , 18.78 , 17.98 , 17.42 , 17.14 , 17.42 , 18.02 , 18.42 , 18.74 , 19.14 , 19.5 , 19])
temporal_window = 42014.0 #seconds

N = len(data) #datapoints
T = 60/(temporal_window/N) #Sample rate average (readings/minute)

###Diagnostic Override###
debug = False #DEBUG SWITCH
if debug:
    wave_lenght = 800 #in seconds (eg. 60*60*2 = 2 hours)
    print "Created a sine wave with %s second period" % wave_lenght
    diagnostic_array = np.arange(0,1,1./N)
    diagnostic_array = np.cos(2*np.pi*temporal_window/wave_lenght*diagnostic_array)
    data = diagnostic_array
#########################
    
a=np.abs(fft.rfft(data, n=data.size))
a[0]=0 #Not sure if this is a good idea but seems to help with choppy data..
freqs = fft.rfftfreq(data.size, d=1./T)
freqs = np.divide(60,freqs)

max_freq = freqs[np.argmax(a)]
print "Peak found at %s second period (%s minutes)" % (max_freq, max_freq/60)

plt.subplot(211)
plt.plot(freqs,a)
plt.subplot(212)
plt.plot(np.linspace(0,temporal_window,data.size),data)
plt.show()

运行上面的代码会产生这个打印语句：

>>>#Using data from database    
Peak found at 1710.49363868 second period (28.5082273113 minutes)

更新

我重写了诊断脚本以进一步测试此代码的可靠性。它允许您创建叠加的波浪集，但也为您提供了一些关于如何表示波浪的选项。

>>>#standing_wave_list = [4,8,9,21,88]
Added a sine wave with 10503.5 second period
Added a sine wave with 5251.75 second period
Added a sine wave with 4668.22222222 second period
Added a sine wave with 2000.66666667 second period
Added a sine wave with 477.431818182 second period
Peak found at 5251.75 second period (87.5291666667 minutes)

如果您想自己尝试，只需剪切并粘贴到之前的代码中即可：

###Diagnostic Override###
debug = True #DEBUG SWITCH
if debug:
    def build_waveform(wave_set, by_period=False):
        #superimposed sine waves (integers create perfet standing waves)
        wave_date = np.zeros(data.size)
        for wave in wave_set:
            if by_period:
                wave_lenght = wave #creates a wave with period n seconds
            else:
                wave_lenght = temporal_window/wave #fits n number of whole waves within the dataset
            new_wave = np.arange(0,1,1./N)
            new_wave = np.cos(2*np.pi*temporal_window/wave_lenght*new_wave)
            wave_date += new_wave
            print "Added a sine wave with %s second period" % wave_lenght
        return wave_date

    option = 2
    if option == 1:
        #####BUILD A SET OF WAVES BY PERIOD IN SECONDS#####
        period_wave_list = [60*60*1,
                     60*30,
                     60*25]
        data = build_waveform(period_wave_list, by_period=True)
        #########
    elif option == 2:    
        #####BUILD A SET OF PERFECT STANDING WAVES#####
        standing_wave_list = [4,8,9,21,88]
        data = build_waveform(standing_wave_list)
        #########
#########################

我保证最终更新！

我发现有必要将 FFT 显示为条形图而不是折线图以使视觉清晰。我还修复了除以零错误（必须在创建时使用“[1:]”语法对数组进行切片）。因此，我将在此处添加代码，但我将删除诊断和数据内容（您可以从以前的代码中复制并粘贴过去）。无论如何，我认为这看起来更清晰：

>>>#standing_wave_list = [4,8,9,21,88]
Added a sine wave with 10503.5 second period
Added a sine wave with 5251.75 second period
Added a sine wave with 4668.22222222 second period
Added a sine wave with 2000.66666667 second period
Added a sine wave with 477.431818182 second period
Peak found at 5251.75 second period (87.5291666667 minutes)

import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt

#data = np.array(just copy and past from previous code)
temporal_window = 42014.0 #seconds

N = len(data) #datapoints
T = 60/(temporal_window/N) #Cycles per minute

###Diagnostic Override###
    #REMOVED
#########################
    
a=np.abs(fft.rfft(data, n=data.size))[1:]
freqs = fft.rfftfreq(data.size, d=1./T)[1:]
freqs = np.divide(60,freqs)

max_freq = freqs[np.argmax(a)]
print "Peak found at %s second period (%s minutes)" % (max_freq, max_freq/60)
plt.subplot(211,axisbg='black')
plt.bar(freqs,a,edgecolor="gray",linewidth=2)
plt.plot(freqs,a, 'r--')
plt.grid(b=True, color='w')

plt.subplot(212,axisbg='black')
plt.plot(np.linspace(0,temporal_window,data.size),data,'r')
plt.grid(b=True,axis="y", color='w')
plt.show()

【讨论】：