6.14. FFT

フーリエ変換は、時間領域の信号を構成周波数成分に分解します。高速フーリエ変換FFT)は、それを効率的に計算するための標準アルゴリズムです。このモジュールが存在する理由であり、また、実用に足る長さの変換(1024点、4096点)がそもそもマイクロコントローラ上で扱える理由でもあります。

カメラ上では、入力は通常、マイク、加速度センサーの軸、電流センサー、または振動プローブから等間隔でサンプリングされたバッファです。出力はスペクトルであり、アプリケーションはそこからピーク、エネルギー、または分類器の特徴を検査します。

numpy.fft は1次元の離散フーリエ変換を提供します。2つの関数は fft()(順変換)と ifft()(逆変換、N で正規化)です。

変換長は2のべき乗でなければなりません。その他の長さは ValueError を発生させます。N = 256N = 1024 がよく選ばれます。ほとんどの組み込みサンプルレートで使える周波数分解能を得るのに十分な長さでありながら、RAMに収まり、1フレーム期間内に十分余裕をもって完了できる短さです。

fft() は、入力の実部を最初の位置引数として、オプションの虚部を2番目の引数として受け取り、実数配列の2要素タプル (real, imag) を返します。

6.14.1. 実数入力、振幅出力

アプリケーションが実数入力のみを与え、振幅スペクトルのみを必要とする場合:

real, imag = np.fft.fft(x)
magnitude  = np.sqrt(real * real + imag * imag)

結果は入力と同じ長さの実数 ndarray です。spectrogram() ヘルパーは、中間的な割り当てなしでまさにこれを行います。ストリーミングループ内ではこちらを呼ぶのが適切です。

6.14.2. 逆変換

ifft() は、往復で浮動小数点誤差の範囲内で元の入力が返るように正規化されています。どちらの呼び出しも実数配列のペアを受け取り、返します:

real, imag       = np.fft.fft(y)
re_back, im_back = np.fft.ifft(real, imag)

6.14.3. 支配的な周波数を見つける

よくあるパターンは、サンプルバッファ(音声、振動、変調IR)の中で最大のピークを見つけることです。手順は次のとおりです:

  1. N を2のべき乗として、N 個のサンプルを取得します。

  2. FFTを実行します。

  3. 最大の振幅を持つビンを見つけます。

  4. サンプルレートを使ってビンのインデックスをHzに変換します。

from ulab import numpy as np

N  = 1024
fs = 8000.0                       # sample rate, Hz

# ... fill ``samples`` with N data points ...

real, imag = np.fft.fft(samples)
spectrum   = np.sqrt(real * real + imag * imag)

half      = spectrum[:N // 2]     # only first half is meaningful
peak_bin  = np.argmax(half)
peak_hz   = peak_bin * fs / N

print("peak at", peak_hz, "Hz")