Kaggle Advent Calendar 22日の記事です。

仕事では宇宙に散らばる銀河の３次元の密度をフーリエ変換して power spectrum を計算するコードの開発とテストをしています。なので、Kaggle にも時々現れるフーリエ変換について基本的なことを少し書いてみます。直接 Kaggle に関係なるわけではないですが、何か「なるほどー」と思うことがあれば幸いです。

図１:（左）銀河の観測データで点ひとつひとつが銀河（右）銀河の３次元の密度分布をフーリエ変換して power spectrum を計算したもの、長さのスケール別の密度ゆらぎの大きさです。密度の高いところは万有引力でどんどん銀河が集まってくるので、ここから宇宙スケールでの重力の法則についてわかります¹。

2019年の Kaggle では

Kaggle でもフーリエ変換が使われることがたまにあります。

電線コンペ (VSB Power Line Fault Detection) では交流電圧データから低周波の交流成分を取り除いて高周波部分のノイズのみに注目するために使われました。
地震コンペ (LANL Earthquake Prediction) では振動の周波数別の振幅である power spectrum が有力な特徴量となりました (時系列では spectral density と言うことが多いかも) 。
音声認識 (Freesound Audio Tagging など) では、短い時間毎に power spectrum を計算して時間と周波数の２次元画像 (spectrogram) にして画像認識するのが定石らしいです。

フーリエ変換とは

$f$ を $t \in [0, T]$ で定義された実数値関数とします。複素数値でもいいですが、実数値データのことが多いと思います。また、時系列データに限らず n 次元空間に対しても同じように使えますが、ここでは１次元の実数データとします。

フーリエ変換とは関数 $f$ から

$$
f(t) = \sum_{n=-\infty}^{\infty} \hat{f}_n e^{i \omega _n t}
$$

のように展開される係数 $\hat{f_n}$ への変換のことをいいます。

i は虚数単位; $i^2 = -1$、
$n$ は整数、
$\omega_n = \frac{2\pi}{T}n$ は角周波数と呼ばれるものです

$i$ の前の符号や定数倍は人や分野によって違うので許してください²。

$e^{i\omega t}$ は

$$ e^{i\omega t} = \cos(\omega t) + i \sin(\omega t) $$
なので三角関数の和で展開していることになります。

逆に係数 $\hat{f}_n$ から f(t) に戻すことを逆フーリエ変換と言います。

図２: 三角関数さん。役に立つことは明らかで三角関数なしでは物理学が始まらないというのに、なぜか時々「社会で役に立たないもの」の代表として弾糾される。

$f(t)$ は実数値関数だったのに $\hat{f_n}$ は複素数、ナニソレコワイ、と思うかもしれません。確かに、

$$
f(t) = \sum_{n=-\infty}^{\infty} A_n \cos(\omega_n t)

B_n \sin(\omega_n t)
$$

のように実数だけで展開してもいいですが、複素数のまま扱うのがお勧めです。

(1) 計算がシンプル、例えば、

$e^{i\alpha} e^{i\beta} = e^{i(\alpha + \beta)} $ ← 簡単！
$\cos(\alpha) \cos(\beta) = $ あー、高校の頃やったよね、$\cos(\alpha + \beta) = \cos \alpha \cos \beta - \sin \alpha \sin \beta$ で... ← cos sin の場合とか sin sin の場合もあって面倒

(2) また、t=0 というのは特別な値ではなくて数値化するときに便宜的上 0 を定めた程度の場合が多いですが、t=0 の位置をずらすと A_n と B_n が混じり合います。つまり、An と Bn は本質的にはひとつのもの（複素数 fn）なのに t=0 の原点の取り方に依存した人為的な分解をされていると言えます。

実数値としてイメージしたい場合は sin cos ではなく

$$ f_n = |f_n| e^{i \phi_n} $$

と極形式で見る方がいいでしょう。t の原点を $\Delta t$ ずらした時,

$$ \hat{f_n} \mapsto \hat{f_n} e^{-i\omega_n \Delta t}$$

と、複素平面で回転します。つまり

絶対値 $|\hat{f_n}|$ は原点 t=0 に依存しない
原点を Δt ずらすと位相が $\phi_n \mapsto \omega_n \Delta t$ ずれる

フーリエ展開の式を極形式で書くと

$$f(t)
= \sum_{n=-\infty}^\infty |\hat{f_n}| e^{i(\omega_n t + \phi_n)}
= \hat{f_0} + 2 \sum_{n=1}^\infty |\hat{f_n}| \cos( \omega_n t + \phi_n )
$$

となります³。信号 f(t) が振幅 2 |f_n| 、角周波数 ωn の cos 波に分解されて、その cos は位相が $\phi_n$ 横にずれているということがわかります。

図３: 位相 φ = -1 の $\cos (\omega t + \phi)$.

離散フーリエ変換

実際には無限にある実数 t すべてでデータは得られないので、等間隔にサンプリングされた
$$ t_m = \frac{m}{N} T \quad （m=0, 1, ..., N-1)$$
で与えられることが多いでしょう。その場合はフーリエ係数は無限にある整数すべての n ではなくて、

$$\hat{f_n}; 0 \le n \le N/2$$

だけあればで十分になります。

$\hat{f_0}$ はいつも実数
N が偶数の時は $\hat{f_{N/2}}$ は純虚数

になるので、N個の実数が N 個の実数に変換されていています。

ではこの範囲外の $\hat{f_n}$ はゼロなのかと言うとゼロではありません。でも、これだけの n があれば必要十分なのです。

(1) まず、n<0 では元の f が実数値ということから、

Reality condition

$$ \hat{f_{-n}} = \hat{f_n}^* $$
（*は複素共役）という式が成り立ちます。なので、n<0 で係数がゼロになるわけではないけれども、独立した情報は持ちません。

(2) フーリエ係数の周期性

$ e^{i \omega_n t} $
と、角周波数 $ \omega_{n + Nl} = 2\pi (n + Nl)/T $ の
$ e^{i [\omega_n + \frac{2\pi N}{T} l ]t} $

というふたつの関数は実数 t 全体で見ると全く異なった関数ですが、等間隔の t = T m/N に限ってい言うと全く同じ値を取ります（l,m は整数）。

図４: N=8 の場合、n=2 の波（青）と n=2 + 8 の波（橙）は等間隔にサンプリングされた時間（灰色の点線）では全く同じ。

なので |n| ≧ N/2 の高周波は離散的な時間では存在しないと言ってよくて、強いてフーリエ係数を計算すると、$e^{i\omega_n t}$ が離散的な t 上で同じなので、同じフーリエ係数が周期的に出てきます。

そんなわけで、離散的有限個のデータではフーリエ係数 $\hat{f_n}$ は N/2 + 1 個の n = 0, ..., N/2 が全ての情報を含んでいます [Nが奇数の場合は (N+1)/2 個]。

高速フーリエ変換 (FFT) と素因数分解

離散フーリエ変換を高速に行うアルゴリズムが Fast Fourier Transform (FFT) です。実数データのフーリエ変換は numpy や scipy の rfft 関数を使ってできます。これで得られたフーリエ係数 $\hat{f}_n$ に逆フーリエ変換 irfft を使うと時間の関数 f(t)に戻ります。

import numpy as np

f = np.ones(1000)      # f(t)
fhat = np.fft.rfft(f)  # フーリエ変換 f(t) -> fhat(omega)
ff = np.fft.irfft(fhat)    # 逆フーリエ変換で f(t) に戻る f == ff
                           # （不動少数点数の範囲で）

基本形はデータの長さが $N=2^n$ の時で計算量は O(N log N) になります。N が２の冪乗でなくても素因数分解して素因数が 3 や 5 だけなら OK です。しかし、あなたが素数大好き東工大生だとしても FFT の N に大きな素数を使うことはおすすめできません。その場合は $O(N^2)$ の計算量がかかってしまいます。timeit で測ってみた FFT の計算時間です。

N	numpy.fft.rfft	scipy.fft.rfft
$1024^2$	30.8 ms ± 443 µs	26.1 ms ± 310 µs
$10^6$	28.3 ms ± 4.35 ms	22.9 ms ± 3.54 ms
$10^6 + 1$	278 ms ± 2.86 ms	194 ms ± 50.3

$10 = 2 \times 5$ は小さい素因数を持っていますが、$N=10^6 + 1$ は良くなくてデータ量はほぼ同じなのに10倍の時間がかかっています。FFTに時間がかかって困った場合はデータの長さもチェックしてみましょう。ところで、Numpy と Scipy でも同じではなく、scipy の方がちょっと最適化されたものを使っているようです。

畳み込み定理 (convolution theorem)

フーリエ変換は Convolutional Neural Network (CNN) でもおなじみの畳み込み(convolution) と相性の良い変換です。畳み込みとフーリエ変換の関係に慣れるともっとフーリエ変換と仲良くなれます。

関数 W と f の畳み込み W*f とは

$$
(W*f)(t) = \int W(t - s) f(s) ds
$$

フィルターだとかカーネルだとか呼ばれる関数 W をスライドしながら掛け算していく演算です。定義域[0, T] は円状になっていて、左端と右側がつながっているとして扱います。

移動平均

W が一部で1、その他は0という単純な関数の時は畳み込みはもっと馴染みのある移動平均になります。

図５: （左）電線コンペのデータ（青）を幅 0.05 秒で移動平均したもの（橙）（右）移動平均に相当するフィルター W(t)。これを移動させながら f(t) と掛け算していく畳み込み積分は移動平均でもあります。

Convolution theorem

$$
\hat{W*f}(\omega) = \hat{W}(\omega) \hat{f}(\omega)
$$

畳み込みした後にフーリエ変換するのは、最初にそれぞれフーリエ変換した後に掛け算したのと同じという定理です。

図６: つまり、電圧データをフーリエ変換した $\hat{f}(\omega)$（左・青）にフィルター関数のフーリエ変換（中央・青）をかけると左・橙になります。これは移動平均したものをフーリエ変換したものと同じで、逆に左図の橙線を逆フーリエ変換すると右図になって移動平均した場合と同じになります。

手間のかかる移動平均の仕方と思うかもしれませんが、FFTは高度に最適化されていて計算が早いので場合によっては畳み込みを計算するよりフーリエ変換を使った方が早い場合があります。

中央の移動平均の W^ を見てみると低周波では 1、高周波では 0 に近く、多少振動はありますが、低周波だけ取り出す low-pass filter の一種だというのがわかります。

Sharp k cutoff

今度は逆にフーリエ係数の高周波成分をサクッと0にして低周波のみ取り出した場合にどうでしょう。

これは $\hat{f}$ に低周波数では1、高周波数では0 になる関数 $\hat{W}(\omega)$ をかけることと同じです。

図７: （左）f^ の高周波成分を0にする（青）ということは、f^ に1か0をとる関数（橙）をかけることと同じ。（右）$\hat{W}(t)$ に対応する W(t).

この W^ を逆フーリエ変換して W(t) を得ると右図の振動する sin(t)/t のような関数になります。つまり、高周波を完全にカットすることはこのフィルター関数で畳み込み積分を行うことと同じになります。このように振動する関数で畳み込みを行うので場合によっては振動のない関数 f(t) に人工的な振動が現れてしまう場合もあります。それが嫌な場合は高周波を突然カットするのではなく、W(t) が振動しない $\hat{W}_n$ で高周波の寄与を減らすのが良いでしょう。ちなみに W(t) が Gaussian のとき、フーリエ変換した W^ も Gaussian になります。時間世界でも周波数世界でも振動しない振る舞いの優雅な関数です。

このように、畳み込み積分をフーリエ係数の世界で見たり、逆にフーリエ係数の周波数別の操作を畳み込み積分で見たりすることでどんな計算がされているのか直感的にわかりやすくなる場合があります。畳み込みのお供にフーリエ変換、フーリエ変換のおつまみに畳み込み積分などいかがでしょうか。

メリークリスマス＆良いお年を。明日は Kaggle Advent Calender の管理者でもあるカレーちゃんさんです。

References

https://arxiv.org/abs/1607.03150, https://arxiv.org/abs/1509.06400 ↩
整合性が取れていればなんでもいいのですが、符号や $2\pi$ などの定数が間違っていたらすみません。ところで、数学科の人に $\pi=3$ などと言えば一生口を聞いてくれなくなること間違いなしですが、「$2\pi = 1$ だよね」ならにっこりしてくれるかもしれません。 ↩
Reality condition と呼ばれる $\omega_{-n} = \omega_n^* $ （* は複素共役）を使いました。フーリエ展開の両辺の複素共役をとることで示せます。 ↩