Python Scipy FFT用于确定具有更高频率信号的间隔

Question

我有一个我正试图解决的问题。我有一个深度，GR，相值，代表取样的深度，伽马值和岩石类型。 我试图设置一个移动窗口（30英尺、15英尺以上和15英尺以下），并找出哪个层段的分层/互层频率最高。互层由伽马射线值反映的不同岩石类型决定

所以我想我会在伽马射线上使用scipy fft来确定每个移动间隔的频率和振幅，然后计算所有频率的p90值（？？我不知道这是不是正确的方法）。我所期望的是一条连续的曲线，指示最高的频率值，但我得到了其他东西。有几个nan（+/-）。这是我到目前为止的代码

from scipy import fftpack
import numpy as np
import pandas as pd
logs_dict ={'Depth': [7075.0, 7075.5, 7076.0, 7076.5, 7077.0, 7077.5, 7078.0, 7078.5, 7079.0, 7079.5, 7080.0, 7080.5, 7081.0, 7081.5, 7082.0, 7082.5, 7083.0, 7083.5, 7084.0, 7084.5, 7085.0, 7085.5, 7086.0, 7086.5, 7087.0, 7087.5, 7088.0, 7088.5, 7089.0, 7089.5, 7090.0, 7090.5, 7091.0, 7091.5, 7092.0, 7092.5, 7093.0, 7093.5, 7094.0, 7094.5, 7095.0, 7095.5, 7096.0, 7096.5, 7097.0, 7097.5, 7098.0, 7098.5, 7099.0, 7099.5, 7100.0, 7100.5, 7101.0, 7101.5, 7102.0, 7102.5, 7103.0, 7103.5, 7104.0, 7104.5, 7105.0, 7105.5, 7106.0, 7106.5, 7107.0, 7107.5, 7108.0, 7108.5, 7109.0, 7109.5, 7110.0, 7110.5, 7111.0, 7111.5, 7112.0, 7112.5, 7113.0, 7113.5, 7114.0, 7114.5, 7115.0, 7115.5, 7116.0, 7116.5, 7117.0, 7117.5, 7118.0, 7118.5, 7119.0, 7119.5, 7120.0, 7120.5, 7121.0, 7121.5, 7122.0, 7122.5, 7123.0, 7123.5, 7124.0, 7124.5, 7125.0, 7125.5, 7126.0, 7126.5, 7127.0, 7127.5, 7128.0, 7128.5, 7129.0, 7129.5, 7130.0, 7130.5, 7131.0, 7131.5, 7132.0, 7132.5, 7133.0, 7133.5, 7134.0, 7134.5, 7135.0, 7135.5, 7136.0, 7136.5, 7137.0, 7137.5, 7138.0, 7138.5, 7139.0, 7139.5, 7140.0, 7140.5, 7141.0, 7141.5, 7142.0, 7142.5, 7143.0, 7143.5, 7144.0, 7144.5, 7145.0, 7145.5, 7146.0, 7146.5, 7147.0, 7147.5, 7148.0, 7148.5, 7149.0, 7149.5, 7150.0, 7150.5, 7151.0, 7151.5, 7152.0, 7152.5, 7153.0, 7153.5, 7154.0, 7154.5, 7155.0, 7155.5, 7156.0, 7156.5, 7157.0, 7157.5, 7158.0, 7158.5, 7159.0, 7159.5, 7160.0, 7160.5, 7161.0, 7161.5, 7162.0, 7162.5, 7163.0, 7163.5, 7164.0, 7164.5, 7165.0, 7165.5, 7166.0, 7166.5, 7167.0, 7167.5, 7168.0, 7168.5, 7169.0, 7169.5, 7170.0, 7170.5, 7171.0, 7171.5, 7172.0, 7172.5, 7173.0, 7173.5, 7174.0, 7174.5, 7175.0], 'GR': [35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 35, 35, 35, 35, 110, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 55, 55, 55, 140, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15, 15, 135, 15], 'Facies': [1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1]}

logs=pd.DataFrame(logs_dict)
target_window = 30
logs=logs.reset_index()
data_freq=logs[['GR']]
for i in range(0, len(logs)):
    if i <= target_window:
        logs.at[i, 'Interbedding'] =0
    elif i > (len(logs)-target_window):
        logs.at[i, 'Interbedding'] = 0
    else:
        data =data_freq.iloc[int(i-target_window/2):int(i+target_window/2),].values
        time_step=data[1]-data[0]
        ps = np.abs(fftpack.fft(data))
        freq = fftpack.fftfreq(data.size,time_step)
        freq=freq[:data.size//2]
        value=np.percentile(freq, 90)*2000
        logs.at[i, 'Interbedding'] =value

facies_colors = ['darkviolet', 'lightgreen', 'r','yellow','navy','coral']
cmap_facies = colors.ListedColormap(
    facies_colors[0:len(facies_colors)], 'indexed')
cluster=np.repeat(np.expand_dims(logs['Facies'].values,1), 100, 1)

#Plotting the results 
f, ax = plt.subplots(nrows=1, ncols=3, figsize=(5, 9))
f.subplots_adjust(top=0.75, wspace=0.1)
for i in range(len(ax) - 1):
    ax[i].set_ylim(top, bottom)
    ax[i].invert_yaxis()
    ax[i].grid()
    ax[i].locator_params(axis='x', nbins=3)
ax00=ax[0].twiny()
ax00.plot(logs.GR, logs.Depth, '-g', alpha=0.8, lw=0.9)
ax41 = ax[1].twiny()
ax41.plot(logs.Interbedding, logs.Depth, '-k', alpha=0.8, lw=0.9)
ax61 = ax[2].twiny()
ax61.imshow(cluster, interpolation='none', aspect='auto',cmap=cmap_facies, vmin=0, vmax=5)
ax00.spines['top'].set_position(('outward', 0))
ax00.set_xlabel('GR', color='green')
ax00.set_xlim(0, 150)
ax41.set_xlabel("Interbedding",color='b')
ax41.set_xlim(0,25)
ax61.spines['top'].set_position(('outward',0))
ax61.set_xlabel('Facies')
ax[1].set_yticklabels([]); ax[2].set_yticklabels([])
ax[1].set_xticklabels([]); ax[2].set_xticklabels([])
ax[0].set_xticklabels([])
plt.show()

这就是我得到的。我所看到的是一条连续的曲线，它会向我指出哪个区间的频率更高。我哪里做错了？我没有信号分析的背景。我想看看我是否能在这里应用这个概念。谢谢你的帮助