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

Question

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

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

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

“时间窗口”，从第一次测量到最后一次测量的时间（秒）。所以在哪里：

T = temporal_window/N #should equal roughly 120 seconds

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

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的峰值随着波长的增加而变小（完全是预期的）。但我不确定的是如何改变它，以便在本例中，峰值与波长匹配（以秒为单位）。有可能用快速傅里叶变换吗？我在读关于IFFT转换回时域的文章，但是对这个主题没有很好的理解，我有点不知所措。。

任何关于如何做到这一点的想法或想法都将非常感谢！！ 如果我没有解释清楚我的意图，请让我知道，我很乐意补充细节。 非常感谢！！

Answer 1

多亏了霍布斯的一点推动，我重新评估了我所看到的东西。

经过进一步的研究，我发现rfftfreq函数与linspace相比非常方便。

所以这是更新后的代码，看起来工作正常。作为一个注释，我得到了“RuntimeWarning:divide by zero”，在这里我做np.divide（60，freques）。然而，这似乎并没有影响结果。

我确实注意到，在脚本的当前诊断部分，它允许在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()