processingの備忘録　-逆離散フーリエ変換-

こんにちは。
今回は「逆離散フーリエ変換」を processing でできるようにしてみようと思います。

これまでに、フーリエ変換シリーズとして何記事か書いてきました。

processingの備忘録 -離散フーリエ変換- - プログラミングの備忘録

processingの備忘録 -高速フーリエ変換- - プログラミングの備忘録

フーリエ変換を利用して絵を描く - プログラミングの備忘録

音声ファイルを扱ったときには、おまけとして音声のフーリエ変換を試してみるということもしてみました。

processingの備忘録 -音声- - プログラミングの備忘録

というわけで、フーリエ変換の逆変換に触れてみようというのが今回の目的です。
中でも、プログラミング向きであろう「逆離散フーリエ変換」というものを processing に実装してみます。

逆離散フーリエ変換とは
実装
実際に変換してみる
まとめ

逆離散フーリエ変換とは

その名の通り、離散フーリエ変換の逆です。

フーリエ変換では時間領域のデータが周波数領域のデータに変換されましたが、逆変換では周波数領域のデータが時間領域のデータに変換できます。
つまり、こんな周波数の波が合わさってできた音です、という情報から実際の音の波形を復元できる、といったところでしょうか。

一般的なフーリエ変換は以下のように表されます。

$$ 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$ は任意の実数で、それぞれ時間と周波数の次元を持ちます。

この逆ということですが、$F(f)$ を入れると $f(t)$ を返してくれるような関数になります。

$$ f(t) = \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} F(f) e ^{ \frac{2 \pi i f t}{N} } $$

これが逆離散フーリエ変換です。
導出に関しては省略しますが、割とあっさり求められるようです。

得られた $f(t)$ は複素数となると予想されますが、フーリエ変換前の $F(f)$ が実数値だった場合、虚部は無視できるということになります。

実装

改めて、逆離散フーリエ変換の式を見てみます。

$$ f(t) = \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} F(f) e ^{ \frac{2 \pi i f t}{N} } $$

指数の部分に虚数単位が入っており、このままでは計算しづらいです。
そこで、あのレオンハルト・オイラーの名を冠した、

$$ e ^{\pm i \theta} = \mathrm{cos} \theta \pm i \mathrm{sin} \theta $$

という数式を利用して、実部と虚部に分けます。

$$ \begin{align} f(t) &= \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} F(f) \left( \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) + i \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \right) \\ &= \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} F(f) \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) + i \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} F(f) \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \end{align} $$

ここで、離散フーリエ変換の式を見てみます。

$$ F(f) = \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) $$

符号や $N$ で割るかが違うくらいで、残りはほとんど同じです。

ということで、離散フーリエ変換のときにつくった関数を割と使えます。
ただし受け取る値 $F(f)$ は複素数なので、実部と虚部それぞれについて計算しなければいけません。

$F(f) = \mathrm{Re} + i \mathrm{Im}$ とおくと、

$$ \begin{align} f(t) &= \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} (\mathrm{Re} + i \mathrm{Im}) \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) + i \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} (\mathrm{Re} + i \mathrm{Im}) \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \\ &= \displaystyle \frac{1}{N} \sum _{f = 0} ^{N - 1} \left( \mathrm{Re} ~ \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) - \mathrm{Im} ~ \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \right) + i \frac{1}{N} \sum _{f = 0} ^{N - 1} \left( \mathrm{Im} ~ \mathrm{cos} \left( \frac{2 \pi f t}{N} \right) - \mathrm{Re} ~ \mathrm{sin} \left( \frac{2 \pi f t}{N} \right) \right) \end{align} $$

と変形できます。

長くなりましたが、これをプログラムにすると以下のようになります。
（inverse discrete Fourier transform の略で IDFT です。）

void idft(){
  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] += (re[j] * cos(2 * PI * f[j] * t[i] / N) - im[j] * sin(2 * PI * f[j] * t[i] / N)) / N;
      imag[i] += (im[j] * cos(2 * PI * f[j] * t[i] / N) + re[j] * sin(2 * PI * f[j] * t[i] / N)) / N;
    }
    
    f[i] = real[i];
  }
}

前述の通り、逆フーリエ変換後の $f(t)$ は実部だけ考えれば良いので f[i] = real[i]; となっています。

実際に変換してみる

これにて逆離散フーリエ変換ができるようになったので、何らかの波をフーリエ変換した後、逆フーリエ変換で元に戻るかを試してみます。

フーリエ変換する波 signal は、40 Hzの正弦波と50 Hzの余弦波の合成波にしてみました。

signal[i] = cos(40 *2*PI*t[i]/N) + sin(50 *2*PI*t[i]/N);

実行自体は簡単で、波の情報を入れておいた signal に対して dft() を実行し、その結果 re と im に対して idft() を実行する、という感じです。

描画に関しては、signal (変換前) と sqrt(re*re+im*im) (フーリエ変換後)、f (逆変換後) を良い感じに図示するようにします。

結果、

ここでは上から順に signal と sqrt(re*re+im*im)、f となっています。
見てみると、一番上と一番下の波形が同じ形になっており、逆変換によって周波数の情報から波形を復元することができました。
(ちなみに、真ん中の図でシグナルが見えている位置を出力するとちゃんと40 Hzと50 Hzになっていました。)

実際の音声処理がどうかは知りませんが、フーリエ変換後のシグナルから特定の周波数を除いたりして逆変換をすると、音声の編集 (?) ができます。
振幅の小さい成分を除くとか、逆に加えることもできると思います。