こんにちは。
今回は「離散フーリエ変換」を processing でできるようにしてみようと思います。
フーリエ変換と聞いてもパッと日常生活における利用例が思いつきませんが、化学系の研究を行っている私からすると、核磁気共鳴 (NMR) 分光や赤外分光 (FT-IR)、X線回折などでよく聞くものであります。
（実際にフーリエ変換をするのは機械なので、普段は変換の結果の方に注目していますが…）
そういえば、たまに何かのMVで音楽に反応して棒グラフが動くみたいな演出があるかと思いますが、おそらくそれもフーリエ変換を利用しています。
ということで、細かい中身は数学科でやってもらうとして、さっくりとですがフーリエ変換に触れてみようというのが今回の目的です。
中でも、プログラミング向きであろう「離散フーリエ変換」というものを processing に実装してみます。
目次
離散フーリエ変換とは
フーリエ変換によって、時間領域のデータが周波数領域のデータに変換できます。
ある関数を三角関数に分解する、という言い方もできると思います。
分光法の観点から言えば、吸収や発光といった現象からその周波数やエネルギーを求めることができる、といった感じでしょうか。
こちらの動画を見ると、フーリエ変換がもっとわかるようになるかもしれません。
【大学数学】フーリエ解析入門①(フーリエ級数展開 I)/全5講【解析学】 - YouTube
では離散フーリエ変換についてですが、その名の通り、離散化したフーリエ変換です。
一般的なフーリエ変換は以下のように表されるようです。
$$ F(f) = \displaystyle \int _{-\infty} ^{\infty} f(t) e ^{-2 \pi i f t} $$
これは連続的で、積分や無限大が入ってきたりと、コンピュータでは扱いづらいものになっています。
そこで、有限の区間だけにして一定間隔ごとに区切っていくことを考えます。
$$ F(f) = \displaystyle \sum _{t = 0} ^{N - 1} f(t) e ^{ -\frac{2 \pi i f t}{N} } $$
ここで、$N$ はデータの個数、$i$ は虚数単位、$t$ や $f$ は任意の実数で、それぞれ時間と周波数の次元を持ちます。
これが離散フーリエ変換です。
先程紹介した動画の1本目を見ていただけたでしょうか。
そこでは、時間領域の関数を原点に"巻き付け"、その重心を考えることで周波数領域の関数に変換できる、と説明されていました。
つまり複素数平面で考えたとき、得られた値の絶対値が周波数領域での大きさに、角度が位相に相当します。
実装
改めて、離散フーリエ変換の式を見てみます。
$$ F(f) = \displaystyle \sum _{t = 0} ^{N - 1} f(t) e ^{ -\frac{2 \pi i f t}{N} } $$
指数の部分に虚数単位が入っており、このままでは計算しづらいです。
そこで、あのレオンハルト・オイラーの名を冠した、
$$ e ^{\pm i \theta} = \mathrm{cos} \theta \pm i \mathrm{sin} \theta $$
という数式を利用して、実部と虚部に分けます。
$$ \begin{align} F(f) &= \displaystyle \sum _{t = 0} ^{N - 1} f(t) \left( \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) - i \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \right) \\ &= \displaystyle \sum _{t = 0} ^{N - 1} f(t) \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) - i \displaystyle \sum _{t = 0} ^{N - 1} f(t) \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \end{align} $$
複素数の型があればこのような式変形は必要ないと思いますが、一旦こうすることにしました。
これをコードにすると、以下のようになります。
（discrete Fourier transform の略で DFT です。個人的には DFT と聞くと density functional theory の密度汎関数理論の方を思い浮かべてしまいますが。）
void dft(){ float[] real = new float[N]; float[] imag = new float[N]; for(int i = 0; i < N; i ++){ real[i] = 0; imag[i] = 0; for(int j = 0; j < N; j ++){ real[i] += signal[j] * cos(2 * PI * freq[i] * time[j] / N); imag[i] += signal[j] * -sin(2 * PI * freq[i] * time[j] / N); } //F[i] = sqrt(real[j]*real[j] + imag[j]*imag[j]); } }
続いて、変換したいデータをつくります。
とりあえず単純な余弦波2つの和を考えることにします。
N は 1000、dt は 1 としました。
void wave(){ for(int i = 0; i < N; i ++){ time[i] = dt * i; signal[i] = cos(100 *2*PI*time[i]/N) + cos(300 *2*PI*time[i]/N); }
後は、良い感じに描画できるようにします。
一番上に変換前の関数 wave()、中・下段にそれぞれ青・赤で変換後の関数の虚部・実部を表示するようにしてみました。
freq が 100 と 300 のところに鋭いピークが見られ、変換前の関数の周波数がわかりました。
ここで、計算に使った全ての周波数成分について図示すると、freq が 900 と 700 のところにもピークが見えます。
これは「エイリアシング」とか「ナイキスト周波数」とよばれるものと関連しているそうで、どうやら、標本点が偶然重なってしまうことによってでてきてしまう偽信号らしいです。
よって、変換後の関数の描画は N の半分まででないと正しくないということになります。
まとめ
以上、離散フーリエ変換の計算を processing でできるようになりました。
今回の方法は数式をそのまま実装したような形ですが、これでは $N$ が大きくなると時間がかかってしまいます。
より速い他の方法として「クーリー・テューキーのアルゴリズム」や「ストックホルムのアルゴリズム」などがあるそうです。
実際の分光装置にはこれらのような「高速フーリエ変換 (FFT)」とよばれるものが実装されていると思いますので、そのうち記事にしてみようと思います。
追記 (2023/08/17)
書きました。
以上
それでは、最後まで読んでいただいてありがとうございました。また次回。
