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

Question

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

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

"data" 是从数据库中获取的数据。通常可以假设它们的间隔为 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 以转换回时域的内容，但由于没有很好地理解该主题，我有点不知所措..

任何有关如何实现这一目标的想法或想法将不胜感激!如果我没有清楚地解释我的意图，请告诉我，我很乐意添加详细信息。非常感谢!!

Answer 1

感谢霍布斯的一点插入，我重新评估了我实际看到的内容。

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

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

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

因此，要真正看到其实际效果(其中峰值 FFT 与输入波形周期匹配)，您所要做的就是更改此行:

-来自-

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

-到-

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