フーリエ変換のざっくり解説とPythonによるFFTの実践

時系列データを扱っていると、しばしばフーリエ変換に出会うことがあるかと思います。
しかし、大学とかで一度勉強していないと、なんのこっちゃわからないですよね。
それに実務でやる場合にはPythonとか何かしらプログラミング言語でコードを書かないといけないため、またハードルが高いです。

そこで今回の記事では、まずフーリエ変換がなんたるかのイメージを説明し、そのあとに実際にPythonコードを書いていきます。

今回の記事を読んでいただければ、ざっくりとではありますがフーリエ変換を理解し、Pythonでフーリエ変換のコードを最低限かけるようになるのではないかと思います。

フーリエ変換は理論をしっかり理解しようと思うと、それなりに勉強しないといけませんが、大枠を押さえつつ、簡単なコードを書くだけでしたらそこまで難しくありません。

ぜひ、トライしてみてください。

目次

1. フーリエ変換はどういう時に使う？

まず、フーリエ変換に出会うのは、だいたい時系列データと扱うときではないかと思います。
データを見ていて、そのデータの中に振動成分、つまり周期性がありそうかどうかを知りたい時によく使います。
（画像分野でも使いますが、今回は時系列データでイメージしていきます。）

時系列の波形を見ただけでは、ただのノイズなのか、それとも本当にある振動成分が強いのかどうかなんてわかりませんよね。
それをしっかり調べる方法の一つがフーリエ変換です。

2. どんなグラフもいろいろな周期の三角関数の重ね合わせで表現できる

では、次にフーリエ変換という考え方の大前提を説明します。
このフーリエ変換を考える上で一番重要なのが、どんなグラフもいろいろな周期の三角関数の重ね合わせで表現することができるというものです。

ちょっと例を見てみましょう。

こちらのグラフはなんだかグジャグジャとしたグラフになっていますね。でもよくみるとただのホワイトノイズではなく、どことなく周期性があるように見えます。

実際このグラフは以下の式で書くことができ、4つの三角関数の和で表現できています。

たった4つの周期関数を重ね合わせただけでも、ちょっと複雑なグラフになるのですから、ものすごくたくさんの周期関数の重ね合わせであれば、どんなグラフも作れそうな気がしますね！

最初に「フーリエ変換は周期性がありそうかどうかを知りたい時に使う」と言いましたが、ここまでのことをまとめるとフーリエ変換は「いろいろな周期の三角関数に分解する」と考えることができます。

この「いろいろな周期」が離散的であれば離散フーリエ変換、連続的であればフーリエ変換となります。
そして、コンピュータで計算する際には、離散フーリエ変換になり、この離散フーリエ変換を高速に動かすアルゴリズムが高速フーリエ変換（FFT : Fast Fourier Transform）になります。

なので、実際にPythonでコードを書く際には、このFFTを書くことになります。
ただ、このFFTは偉大な先人がライブラリを作ってくれているので、私たちはほんの数行書くだけで使うことができます。
ありがたやありがたや。

3. 用語の確認

周期・周波数・角周波数

フーリエ変換をやる前に、確認したおきたいワードがあります。
「周期・周波数・角周波数」です。

こちらのグラフは周期＝2π（~6.28秒）です。1つの波が元に戻るまでの時間を周期といいます。

次に周波数は「1秒間に振動する回数」のことなので、6.28秒で1つの波があるのですから、周波数fはf = 1/6.28 [Hz]です。単位はヘルツ：Hzになります。

最後に角周波数ですが、これは「1秒あたりにどれくらいの角度(rad)進むか」なので、角周波数ωはω = 2πf = 1となります。

低周波・高周波

低周波・高周波という言葉がありますが、これは周波数が高いか低いかです。
周波数が高い場合は、「1秒間に振動する回数が多い」という意味です。
一方で周波数が低い場合は「1秒間に振動する回数が少ない」となります。

低周波のグラフの例

高周波のグラフの例

基本的な数学（複素数）

最後の確認として複素数を見ておきましょう。
複素数は実部と虚部を持ち、$${a+bi}$$の形で表すことができます。
これを複素数平面にプロットすると以下のようになります。
複素数平面では、x軸が実部でy軸が虚部になるので、複素数は、$${a+bi=r(\cos⁡ θ+i \sin ⁡θ)}$$と表せます。

また、この複素数は指数で表現することもできます。
$${\cos⁡ θ+i \sin ⁡θ = e^{i\theta}}$$ 
この式はオイラーの公式と呼ばれ、複素数の式に適用すると
$${a + bi =r(\cos⁡ θ+i \sin ⁡θ) = re^{i\theta}}$$ 
になります。

4. フーリエ級数展開

フーリエ変換に入っていく前の導入として、フーリエ級数展開を学んでおく必要があります。
級数展開なので、これはΣを使った周期波形の重ね合わせになるのですが、フーリエ級数展開は周期的な波形を表す理論です。
フーリエ変換には周期性の条件は不要ですが、フーリエ級数展開はあくまで周期波形を、いろいろな周期の関数で表すものであることに注意してください。

フーリエ級数展開の式は以下の通りです。

$${x(t)=\frac{a_0}{2}+\sum_{j=1}^{\infty} a_j \cos(2\pi j \frac{t}{T})+\sum_{j=1}^{\infty} b_j \sin(2\pi j \frac{t}{T})}$$

$${a_j=\frac{2}{T}\int_{t_0}^{t_0+T} x(t) \cos(2\pi j \frac{t}{T})dt}$$
$${b_j=\frac{2}{T}\int_{t_0}^{t_0+T} x(t) \sin(2\pi j \frac{t}{T})dt}$$

この式のなかでaとbという文字がありますが、これらは各周期成分の振幅を表しており、フーリエ級数展開において求めたい変数になります。このaとbは各周波数の波形に対して強弱をつけるイメージです。
Σのjが1から∞になっており、∞ではなくある程度で打ち切るようにすれば、対象の関数を近似するイメージになります。

さらにこのフーリエ級数展開は、複素数を用いて複素フーリエ級数展開として表すこともできます。

$${x(t)=\sum_{j=-\infty}^{\infty} c_j e^{2\pi j \frac{t}{T}i}}$$

$${c_j=\frac{1}{T}\int_{-T/2}^{T/2} x(t) e^{-2\pi j \frac{t}{T}i}dt}$$

$${a_j = 2Re(c_j)}$$
$${b_j = -2Im(c_j)}$$
$${c_j=\frac{a_j-ib_j}{2}}$$

5. フーリエ変換

さて、ここからフーリエ変換をやっていくわけですが、先ほどやったフーリエ級数展開とどう違うのかと思う方も多いと思います。

それは、フーリエ級数展開は周期関数を対象としていましたが、フーリエ変換では非周期関数に対しても扱うことができるのが違いです。

複素フーリエ級数展開では周期関数を対象としており周波数はとびとびの値を用いていましたが、フーリエ変換が非周期関数を対象に扱う場合は、周期Tを∞にすればいいのでは、と考えて式変形をします。

すると、とびとびだった周波数fが連続値になることもわかるかと思います。

このようにフーリエ級数展開が周期関数を対象にしていたことを出発点として、フーリエ変換では周期Tを∞にすると考えることで、フーリエ変換の世界を考えていくことができます。

では、フーリエ変換の式はどうなるかというと、時系列データx(t)に対して、フーリエ変換後の値をX(ω)とすれば

$${X(\omega)=\int_{-\infty}^{\infty} x(t) e^{-i \omega t}dt}$$

のように書くことができます。
この変換はフーリエ変換で、逆にX→xを行う操作をフーリエ逆変換といいます。

では、このように周期Tを無限大にするとほかにどのような変化があるでしょうか？

まず、振幅のスペクトルを見ていくと、フーリエ級数展開では周期がとびとびでしたので、以下のようなイメージです。

しかし、フーリエ変換になると周波数fが連続値になるわけですから、イメージとしては以下のようなスペクトルになります。

6. 離散フーリエ変換

フーリエ変換にて非周期関数を扱うことができるようになったわけですが、まだまだコンピュータで扱うには問題があります。
それは、結局データはとびとびの値であり、コンピュータもとびとびの値しか扱えないというところです。でも非周期関数を扱いたい。
なので、実際にコンピュータで扱うのは、フーリエ変換を離散化した離散フーリエ変換になるということです。

式は以下のようになります。

$${X(k)=\sum_{n=0}^{N-1} x(n) e^{ \frac{-2\pi kn}{N}i}}$$

一旦、フーリエ級数展開とフーリエ変換、離散フーリエ変換の特徴をまとめると以下のようになります。

そしてこの離散フーリエ変換を高速で行うアルゴリズムがFast Fourier Transform (FFT)です。

Pythonやいろいろなプログラミング言語でライブラリが開発されているので、簡単に使うことができます。