TL;DR

サンプル時間が一定でない信号のフーリエ変換を簡易的に行う
- scipyのinterp1dで固定サンプルの信号へと変換
- あとは通常通りnumpyのfft
ちゃんとした解析をしたい場合はPyNuFFTとかを使いましょう（きちんとした理論があります。）

依存パッケージ

numpy, scipy, matplotlibを使います。

from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

0. Test用データ作成

既にデータがある人はSkipしてOK。
ここではrandomとcumsumを使って適当なランダムサンプリング時間の信号を作成します。

# randomなサンプル時刻を作成 200点数
dx = np.random.rand(200)/20 # 0~1の乱数/20≒期待値的には40Hzくらいの信号
x = np.cumsum(dx)

# 正弦波信号とする
f = 3. # hz
y = np.sin(2*np.pi*f*x)

# Plot
plt.plot(x,y)

ちょっと気持ち歪んでるかな？くらいの信号になります。

これのサンプル時間を可視化すると下記のような感じ。

1. scipy.interpolateを使って固定サンプルの信号に変換

x軸を固定サンプル時間にしてそこのyにあたる量を元の信号から補間をします。
下記のコードではわかりやすさのためにサンプル時刻や時刻の範囲に固定値を使ってますが、本当は元の信号から抜き出して来たほうが良いです。

# 新しい時系列の設定
new_sampling_freq = 30. # resampling後のFreq
new_dt = 1./new_sampling_freq
time = np.arange(0,4,new_dt) # 0s-4sまでの固定サンプル時刻


# 補間
func = interpolate.interp1d(x, y, "cubic", fill_value="extrapolate") # f(x)を作成, cubic補間

new_y = func(time)

# 
plt.plot(time, new_y)

サンプリングが揃うので見た目もちょっと揃います。

2. （numpyで）FFT

scipyでFFTしても良いですが、numpy版に良さげな記事があったのでそれに準拠しておきます。

# fft
N = len(new_y)
fft_freq = np.fft.fftfreq(N, d=new_dt)
y_fft = np.fft.fft(new_y)
fft_amp = abs(y_fft/(N/2))

# plot
plt.plot(fft_freq[1:int(N/2)], fft_amp[1:int(N/2)]) # A-f グラフのプロット
plt.xscale("log") # 横軸を対数軸にセット
plt.grid()

もともと設定していた3Hzの信号がしっかり検出できているのがわかります。

関数化まとめ

ということで一連の流れを関数化するとこんな感じになります。


"""
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt
"""


def approx_nufft(x, y):
  if not len(x) == len(y):
    print("length of x and y must be the same.")
    return 
  N_ = len(x)
  max_x, min_x = max(x), min(x)
  mean_dx = (max_x - min_x) /N_
  
  # 1. interpolatoin
  x_new = np.arange(min_x, max_x, mean_dx)
  f_new = interpolate.interp1d(x, y, "cubic", fill_value="extrapolate")
  y_new = f_new(x_new)

  # 2. FFT
  N = len(x_new)
  fft_freq = np.fft.fftfreq(N, d=mean_dx)
  y_fft = np.fft.fft(y_new)
  fft_amp = abs(y_fft/(N/2))

  # plot
  plt.plot(fft_freq[1:int(N/2)], fft_amp[1:int(N/2)]) # A-f グラフのプロット
  plt.xscale("log") # 横軸を対数軸にセット
  plt.grid()
  plt.xlabel("Freq [Hz]")

# test
approx_nufft(x,y)

補足（？）

上記関数では自動でリサンプリングの粗さを設定していますがサンプルの粗さによっては下記のようにちょっとFFTの結果が理想的でなくなる場合もあります。

サンプル時間が一定でない信号のフーリエ変換を近似的に得る(Python)