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);
  }
}