はじめに
Pythonには高速フーリエ変換が簡単にできる「FFT」というパッケージが存在します。
とても簡便な反面、初めて扱う際にはいくつか分かりにくい点や注意が必要な点がありました。
ということで、私がつまづいた箇所などを踏まえて、FFTの使い方をできるだけ分かりやすく本記事にまとめておこうと思います。（クドいほど強調した箇所もあるので、必要がなければ読み飛ばしてもらって構いません。）
用語
FFTとは
FFTとは、高速フーリエ変換のことです。
しかし、Pythonの FFT を扱うにあたって、高速フーリエ変換の中身を学ぶ必要はありません。
実装レベルでは、「離散時間フーリエ変換」「離散フーリエ変換」までの知識があれば素晴らしいですが、そこまで必要ではないと思います。「フーリエ級数展開」「フーリエ変換」はもちろん押さえておいてください。
ちなみに高速フーリエ変換として知っておくべきことは、次のことです。
FFT（高速フーリエ変換）
サンプル数 $N$ が2の累乗の場合のみ、高速フーリエ変換される。
そうでない場合は、離散フーリエ変換が行われる。
高速フーリエ変換を使うと、離散フーリエ変換では $O(N^2)$ の時間かかってしまう計算が、$N = 2^k$ の場合にのみ、$O(N\log N)$ の時間で計算できます。
今回は $N$ がそこまで大きくないので、自由に $N$ を指定しています。
また、以下にも注意してください。
高速フーリエ変換は、離散フーリエ変換の仲間である。
(高速フーリエ変換) \subset (離散フーリエ変換)
離散フーリエ変換と全く別の理論、というわけではありません。（だいぶアクロバティックで異次元なことはしてますが…）
「離散フーリエが分かれば大丈夫」と言うわけはここにあります。
サンプリングレート（サンプリング周波数）
サンプリングレートとは、「波形をデジタルで扱う際に、1秒間あたりいくつサンプリングするか」を表した値のことです。
例
- CDのサンプリングレートは基本的に 44.1kHz である。
- 最近のストリーミング音楽では 48kHz のものも増えてきている。
今回はサンプリングレート $f_s = 44.1$ [kHz] とします。
FFTの使い方
そのまえに…
高速フーリエ変換は、$N$個のデータを、別の$N$個のデータに変換するものです。
これからの説明では、便宜上 $y(t)$ のように関数として説明しますが、実際は $y_i \; (i=1, 2,..., N)$ のように離散化されています。そこだけ注意して読んでいただけると幸いです。
波のグラフ y(t)
では本題に参りましょう。
まず、時間に関する波の式 $y(t)$ を準備します。今回は簡単な例として振幅が1のサイン波
y(t) = \sin(2\pi f_0 t)
を用います。これを音波の式と思うことにすると、$f_0$ は（一定の）周波数です。
このグラフをPythonくんにお願いして描いてもらっちゃいます！
Python初心者のためにゴリゴリにメモをしておきました🦍📜
import numpy as np # numpyをインポート
import matplotlib.pyplot as plt # matplotlibをインポート
### 量子化(離散化) ###
f_s = 44100 # サンプリングレート f_s[Hz] (任意)
t_fin = 1 # 収録終了時刻 [s] (任意)
dt = 1/f_s # サンプリング周期 dt[s]
N = int(f_s * t_fin) # サンプル数 [個]
### 入力信号 y(t) ###
f0 = 440 # 周波数 f0[Hz]
t = np.arange(0, t_fin, dt) # 離散時間 t[s]
y = np.sin(2*np.pi*f0*t) # y(t) を生成
### y(t) のグラフ ###
plt.plot(t, y) # y-t グラフのプロット
plt.xlim(0, 10*10**-3) # 横軸に関する描画範囲指定
plt.show() # グラフの表示
これを実行すると、以下のグラフが得られます。
これが周波数 440 Hz、ドレミで言うと「ラの音」の波形です🎶
補足
サンプリング周期 $dt$
サンプリング周期 $dt$ [s] は、サンプリング周波数 $f_s$ を用いて、
dt = \frac{1}{f_s}
と表せる。
サンプル数 $N$
サンプル数 $N$ は、$\newcommand{\tfin}{t_\mathrm{fin}} \tfin$ と $dt$ を用いて、
N = \frac{\tfin}{dt} = f_s \tfin
と表せる。
波のグラフを周波数のグラフに（フーリエ変換）
音の正体を知らずにこの波形だけを見たとき、「コイツは本当に 440 Hz の音なのか？」と疑いたくなりますよね。
そこで登場するのが FFT です。
FFTを実行することで、音のスペクトルのグラフを得ることができます！
音楽をやっている人なら、"イコライザー(EQ) のグラフ" として見たことがあるかもしれません。
（あるいは、"Photoshopのあのグラフ" として見たことが以下略…）
FFTの実行と、スペクトルを得るプログラムは以下のようになります。
### FFT: tの関数をfの関数にする ###
y_fft = np.fft.fft(y) # 離散フーリエ変換
freq = np.fft.fftfreq(N, d=dt) # 周波数を割り当てる（※後述）
Amp = abs(y_fft/(N/2)) # 音の大きさ（振幅の大きさ）
やっていることは単純で、
- 時間 $t$ の関数 $y(t)$ をフーリエ変換し、周波数 $f$ の関数 $\newcommand{\ty}{\tilde y_\mathrm{fft}} \ty(f)$ を得る。
- サンプル数 $N$ と、$f$ を関係づける。$\ty(f)$ の周波数 $f$ には
freqに格納した値を用いる。 - 各周波数 $f$ ごとの音の大きさ $A$ を算出する。
の3つです。この3ステップを詳しく見ていきましょう。
1. フーリエ変換
フーリエ変換は
\ty(f) = \mathcal{F}[y(t)]
と書けるのでした。（$\mathcal{F}$ は、「 $[\quad ]$ 内の関数をフーリエ変換しますよ～」という演算子です。）
具体的な表式（$e^{-i\omega t}$とか）は今回あまり意識しなくてもOKです。
Python では、 np.fft.fft(y) のように書き、y 部分に変換したいデータを入力します。
とにかく、時間の関数が周波数の関数に変換されることが大事です。
2. 周波数 f の関係づけ
離散フーリエ変換を学ぶと分かりますが、フーリエ変換で生まれた関数 $\ty(f)$ の 周波数 $f$ は、実は $t$ と同様に離散化されています。これにともなって、$f$ は $N, dt$ に依存しているのです。
つまり、$f$ は $N, dt$ の値（すなわちサンプル数が多いか少ないか）によって関係が変化します。固定された $N$ の場合、$f$ は、サンプルの $i$ 番目の関数$(= f(i))$ です。
サンプルの数 $N$ とサンプリング周期 $dt$ によって、 $f$ が変わってくる。
なので、np.fft.fftfreq(N, d=dt) にサンプル数 N とサンプリング周期 dt を代入する、というわけなんですね。
また、$f/2$ [Hz] の周波数を、ナイキスト周波数 と呼び、それより大きい周波数 $f$ に関して、$\ty(f)$ はオバケのような値をとります。理論的に変な値になってしまうことが知られているということです。したがって、グラフで表す際はその値を除いてやる必要があります。
今回、サンプリングレートを $f_s = 44100$ [Hz] に設定したので、 $f = 22050$ [Hz] がナイキスト周波数となります。
ナイキスト周波数： $f_s/2$
これより大きい周波数では怪奇現象が起こる👻（や、ちゃんとした理論があります）
3. 音（振幅）の大きさ
PythonライブラリnumpyのFFTの仕様では、$f = f_1$ 地点での振幅 $A$ は、次のように表されます。
$f = f_1$ 地点での振幅 $A$
A = \left| \frac{\ty(f_1)}{N/2} \right|
ナイキスト周波数以下のデータ数 $N/2$ で割ることにより、$y(t)$ の振幅と一致します。
また、$\ty(f_1)$ は複素数なので、絶対値をつけなければなりません。
音波のスペクトル
前節《周波数 f の関係づけ》でも説明したように、$\ty(f)$ のうち、意味のあるデータというのは、ナイキスト周波数以下のデータのみです。つまり、各データの（1から数えて）$2$ 番目から $N/2$ 番目までのデータが意味のある値です。
したがって、Pythonの基本技「スライス」を利用して、グラフに用いるデータを制限します！
Pythonでの処理
周波数のデータ freq を、（0から数えて） 1 番目から int(N/2) - 1 番目までのデータにスライスすると、 freq[1:int(N/2)] となり、振幅のデータ Amp に対しても同様に、 Amp[1:int(N/2)] となります。
ここで、 0 番目のデータを除く理由は、次のとおりです。
0 番目のデータは 0 Hz なので、使わない！
以上の点に注意して、スペクトルのグラフを描くと以下のようになります。
### 音波のスペクトル ###
plt.plot(freq[1:int(N/2)], Amp[1:int(N/2)]) # A-f グラフのプロット
plt.xscale("log") # 横軸を対数軸にセット
plt.show()
これを実行すると、
のグラフが得られますが、分かりにくいので拡大すると、
のようになっています。ちゃんと $f = 440$ [Hz] のところに、振幅 $A = 1$ のピークが立っていますね。
フーリエ逆変換でもとの波形に戻す
$\newcommand{\iy}{y_\mathrm{ifft}}$フーリエ逆変換は、
\iy(t) = \mathcal{F}^{-1}[\tilde y(f)]
で表せました。これも具体的な表式（$e^{i\omega t}$とか）を意識しなくて大丈夫です。
Python では、 np.fft.ifft(y) のように書いて、y 部分に変換したいデータを入力します。
今回は、$\ty$ を $\iy$ に変換したいから、
### IFFT: fの関数をtの関数にする ###
y_ifft = np.fft.ifft(y_fft)
と入力すればよいですね。
このグラフを出力するには、以下のコードを食わせましょう。
### 逆変換したグラフ: もとの波形を復元 ###
plt.plot(t, y_ifft.real) # 実部しかいらない
plt.xlim(0, 10*10**-3)
plt.show()
このとき、$\iy$ は実部をとって、 y_ifft.real とすることに注意。（ダッテ私タチ実数世界ニ住ンデルカラ…虚数ノ波ッテ見エナイデショ？）
これを実行すると、
となり、$\iy(t)$ が、最初に設定した $y(t)$ と一致したことがわかります。
あらためて
高速フーリエ変換は、$N$個のデータを、別の$N$個のデータに変換するものです。
あたかも滑らかな関数 $y(t)$ かのように説明してしまいましたが、実際は $y_i \; (i=1, 2,..., N)$ のように離散化されています。そこだけ注意して理解してください。
あとがき
これでPythonでのフーリエ変換/逆変換ができるようになりました。お疲れさまでした。
余裕があれば、応用として $y(t)$ が三角波・のこぎり波・矩形波の場合もやってみると面白いと思います。
なにか誤りや質問等ありましたら、気軽にコメントいただけると嬉しいです。
参考文献
コグニカルは理論を知るのに分かりやすいのでおすすめです。