こんにちは。
今回は「逆高速フーリエ変換」を processing で実装してみます。
なかなか長編になってきました、フーリエ変換シリーズです。
processingの備忘録 -離散フーリエ変換- - プログラミングの備忘録
processingの備忘録 -逆離散フーリエ変換- - プログラミングの備忘録
processingの備忘録 -高速フーリエ変換- - プログラミングの備忘録
音声ファイルを扱ったときには、おまけとして音声のフーリエ変換を試してみるということもしてみました。
processingの備忘録 -音声- - プログラミングの備忘録
先日「逆離散フーリエ変換」を実装し、締めの部分で「逆変換の高速化はまだ先になりそうです」と書きました。
が、意外と簡単にできそうだったので「逆高速フーリエ変換」も扱ってみます。
(これに伴い、順変換の記事も理解が不十分だった部分を若干補充しました。)
目次
理論
高速フーリエ変換は過去にやりましたので、その逆も同じように考えていけばいいはずです。
数式自体は似ているので、工夫すれば一から考えなくても良いかもしれません。
しかし、(こう言っては失礼ですが) 面倒です。
一度、フーリエ変換の式を見てみましょう。
離散フーリエ変換
$$ F(f) = \displaystyle \sum _{t = 0} ^{N - 1} f(t) e ^{ -\frac{2 \pi i f t}{N} } $$
逆変換
$$ f(t) = \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} F(f) e ^{ \frac{2 \pi i f t}{N} } $$
逆変換の式で右辺全体の複素共役をとってみます。
同時に、等しくなるように細かな部分をいじります。
$$ f(t) = \displaystyle \frac{1}{N} \overline{ \sum _{f = 0} ^{N - 1} \overline{F(f)} e ^{ -\frac{2 \pi i f t}{N} } } \\ $$
すると、$F(f)$ の複素共役のフーリエ変換の複素共役を $N$ で割ったものが逆変換に等しい、ということがわかりました。
つまり、わざわざ逆変換を考える必要はなく、既にあるフーリエ変換の関数を使い回すことができるのです。
(逆離散フーリエ変換でも同じ考え方ができるはずですが、記事更新当時は知らなかったので数式をそのまま実装していました。)
普通に Wikipedia に書いてありました。
高速フーリエ変換 - Wikipedia
実装
というわけで、実装してみます。
Complex[] IFFT(Complex[] input){ int N = input.length; Complex[] output = new Complex[N]; // F(f) の複素共役 for(int i = 0; i < N; i ++){ input[i] = input[i].conjugate(); } // 高速フーリエ変換 output = FFT(input); // 変換後の値の複素共役を N で割る for(int i = 0; i < N; i ++){ output[i] = (output[i].conjugate()).cdiv_scalar(N); } return output; // f(t) の出力 }
ここでは、複素数のクラスを使っています。
(高速フーリエ変換のときとは少し違っているので気をつけてください。)
class Complex{ float re, im; Complex(float _re, float _im){ re = _re; im = _im; } Complex cadd(Complex a){ return new Complex(re+a.re, im+a.im); } Complex csub(Complex a){ return new Complex(re-a.re, im-a.im); } Complex cmult(Complex a){ return new Complex(re*a.re-im*a.im, re*a.im+im*a.re); } Complex cmult_scalar(float a){ return new Complex(re*a, im*a); } Complex cdiv(Complex a){ float b = a.re*a.re + a.im*a.im; return new Complex((re*a.re-im*a.im)/b, (-re*a.im+im*a.re)/b); } Complex cdiv_scalar(float a){ return new Complex(re/a, im/a); } float cabs(){ return sqrt(re*re + im*im); } Complex cexp(){ return new Complex(exp(re)*cos(im), exp(re)*sin(im)); } Complex conjugate(){ return new Complex(re, -im); } }
試用
20 Hzの正弦波と70 Hzの余弦波の合成波で試してみました。
上から、変換前、変換後、逆変換後、となっています。
変換後は絶対値と実部、虚部をそれぞれ描画しています。
ピークのある位置は20 Hzと70 Hzで、合っているようです。
これを逆変換して、実部と虚部をそれぞれ描画しています。
実部 (赤) は元の形に戻っており、虚部 (青) はほとんどないことがわかります。
どうやら合っていそうな描画となっており、これにて逆高速フーリエ変換の実装完了です。
ちなみに、計算速度は離散化しただけのものよりもかなり速くなっていました。
まとめ
以上で、逆高速フーリエ変換を processing で実装することができました。
高速にできたので、そのうち二次元に拡張して画像処理に使ったり、スペクトログラムをつくってみたりという記事も書くかと思います。
お楽しみに。
最後まで読んでいただいてありがとうございました。また次回。