プログラムdeタマゴ

nodamushiの著作物は、文章、画像、プログラムにかかわらず全てUnlicenseです

高速フーリエ変換(FFT)の解説。理論編

私はあまり画像を波長空間でフィルタリングとか言うことをやらないので、実のところ、今までFFTどころか離散フーリエ変換(DFT)すらしたこと無かった。

というわけで、ちょっと調べてみたのでまとめてみようと思う。かなり長い記事になるよ。

[原理解説]

 離散フーリエ変換の式(1)について考察をおこなう。
F(u) = \sum_{x=0}^{N-1}f(x)e^{\frac{2\pi i}{N}xu}  (1)
 式(1)において、特に
N=A^n \\0\le u \lt N-1\\A,n,u \in Natural
となる場合を考える。なお、nに関しては0を含んでも良い。(Nが全く見た目が変化しなかったので、英語で書いた。要するに自然数)


 このとき、
0 \le j < A^{n-1} \\0\le k < A\\ 0\le l \lt A\\ 0 \le m \lt A^{n-1}
を満たす適当な整数j,k,l,mを定義すると、u,xは以下の様にかける

u=jA+k  (2)
x=lA^{n-1}+m  (3)


(2),(3)を用いて式(1)を書き換えると、
F(jA+k)=\sum_{m=0}^{A^{n-1}-1}\sum_{l=0}^{A-1}f(lA^{n-1}+m) exp(\frac{-2\pi i}{A^n}(jA+k)(lA^{n-1}+m))  (4)
となる。

式(4)を展開し、整理すると、
F(jA+k)=  \sum_{m=0}^{A^{n-1}-1}exp(\frac{-2\pi i}{A^{n-1}}jm)\{ exp(\frac{-2\pi i}{A^n}km) \sum_{l=0}^{A-1}f(lA^{n-1}+m)exp(-\pi ikl)\}  (5)


ここで、式(5)の大括弧でくくった部分は、mとkの関数と見なせるので
\cos\theta+i\sin\theta= W_\theta  (6)
を用いると、結局式(1)は
F(jA+k) = \sum_{m=0}^{A^{n-1}-1}f_1(m,k)W_{Ajm\theta}  (7)
f_1(m,k)=W_{km\theta} \sum_{l=0}^{A-1}(-1)^{kl}f(lA^{n-1}+m)  (8)
\theta =-\frac{2\pi}{A^n}  (9)
という、kを定数と見なせばA-1次のDFTに変換できる。


式(1)を用いてuの定義域の全ての値を求めるのに必要なループ回数はN^2であるのに対し、式(7)を用いて求めた場合はN^2/Aとなる。さらに式(7)を同様に展開していくことで、さらにループ回数を減らすことができる。最終的にはO(NlogN)のオーダーになる。これがFFTの論理背景である。




実用上は式(8)がA=2またはA=4の時単純になるため、どちらかが使われることが多い。ここではA=2の場合として話を進める。



式(8)にA=2を代入すると
f_1(m,k) = W_{mk\theta}\{f(m)+(-1)^kf(2^{n-1}+m)\}

さらに、kは0または1なので、それぞれ計算すると、
f_1(m,0)= f(m)+f(2^{n-1}+m)  (10)
f_1(m,1)= W_{m\theta}\{f(m)-f(2^{n-1}+m)\}  (11)


ここで、f(x)を実部と虚部に分けて
f(x)=R_x+iI_x
と表記する事を考える。R_x,I_xは実数である。


なお、2^{n-1}+mと表記すると長く、読みづらいため、以下の様に略記する。
2^{n-1}+m = m'  (12)


これらを用いて、(10)(11)を表現すると、
f_1(m,0)= (R_m+R_{m'})+i(I_m+I_{m'})  (13)
f_1(m,1)=W_{m\theta}\{R_m-R_{m'}+i(I_m-I_{m'})\}  (14)


式(14)において、
R_m' = R_m-R_{m'}  (15)
I_m' = I_m-I_{m'}  (16)
とし、Wを三角関数表現に直し、整理すると
f_1(m,1)=(R_m'\cos(m\theta) - I_m'\sin(m\theta)) + i(R_m'\sin(m\theta)+I_m'\cos(m\theta))  (17)



(13)(17)により、f_1(m,k)の虚数表現が得られた。ここまでの議論で計算機上で演算を行うために必要な式は出揃った。ここまでの議論は再帰関数等を用いれば実装するのは容易い。

[理論を実装]

特に深く考えずに、上の議論をそのままJAVAソースコードに変換してみる。
realが実数を表す値で、imaginaryが虚数を表す値である。結果はそのまま渡された配列に上書きすることにする。
余計な処理を書かないために、realとimaginaryの長さは同じとし、2のn乗となる数であることとする。実際にはそれを確かめる処理を入れるべきである。

ソースの(9)などは上の論理の話の式番号に当たる。


ソース1

public static void fft(double[] real,double[] imaginary){
  int n = real.length;
  double theta = -2*Math.PI/n;//(9)
  _fft(n,theta,real,imaginary);
}

private static void _fft(
    int n,double theta,
    double[] real,double[] imaginary){
  if(n<=1)return;//n=1の時、(1)の右辺はf(0)であるから、F(0)=f(0)で変わらないので処理なし。
  int nh = n/2;
  double[] ro = new double[nh],io=new double[nh],//F(j*2+1)を格納する
    re = new double[nh],ie = new double[nh];//F(j*2+0)を格納する
  for(int m=0;m<nh;m++){
    re[m] = real[m]+real[nh+m];//(13)
    ie[m] = imaginary[m] + imaginary[nh+m];//(13)
    double Rm = real[m] - real[nh+m];//(15)
    double Im =  imaginary[m] - imaginary[nh+m];//(16)
    double cos = Math.cos(m*theta);//(17)のcos(mθ)
    double sin = Math.sin(m*theta);//(17)のsin(mθ)
    ro[m] = Rm*cos-Im*sin;//(17)
    io[m] = Rm*sin+Im*cos;//(17)
  }
  //ここまででf(m,0)(← re+i*ie)とf(m,1)(← ro+i*io)計算終了
  //次にそれぞれをフーリエ変換すればいいので、_fftに投げる
  _fft(nh,2*theta,re,ie);//(7)
  _fft(nh,2*theta,ro,io);//(7)
  //計算結果をreal とimaginaryに書き出す
  //偶数奇数を順に書き出せばいいので、
  for(int j=0;j<nh;j++){
    real[2*j] = re[j];
    imaginary[2*j] = ie[j];
    real[2*j+1] = ro[j];
    imaginary[2*j+1]=io[j];
  }
}

これを209,7152個のデータを突っ込んで速度を測ったところ、約2秒かかった。ということは、ここではlogの底を2とすると、log2097152=21だから、DDF(O(N^2)のオーダー)だったら、2*2097152/21 =2.3(日)。軽く2日以上はかかる計算になるのか。うわ〜、すげ。
なんだかもう、これで十分な気がしてきたけど、もうちょっと改善してみよう。


まず、上のソースでは、値を書き込み、次に渡す新たな配列を二つ作っているが、よくよくこのソースを見てみると、realやimaginaryの値がその後使われることはない。ということは、いちいち配列を作るなんてメモリの無駄をせずに、はじめから書き込んでいったらいいと言うことになる。


では、F(2j)の値をf(j)が示す配列の場所に、F(2j+1)の値をf(j+2^(n-1))が示す場所(つまり、real[j+nh]とか)に突っ込んでみよう。
そうすると、_fft(nh,2*theta,re,ie)という、F(2j)を計算する部分は、realとimaginaryの前半部分、F(2j+1)は後半部分を渡せばいい。
さらに、その場合気がつくのは、同じ呼び出し深さのメソッドが値を読み込み、書き込む位置は絶対にかぶらない。と、言うことは毎回並び替えして返さなくても、一番最後に並び替えてしまっても問題ないという事じゃないだろうか。


では、並び替えるにはいったいどうすればいいのか?F(t)がいったいどこに入るのだろうかを考える。
tがもし、2^n +1と表せたら(つまり奇数だったら。さらに言えば0ビット目が1だったら=>t&1が1だったら)配列の後半に来るはずだ。逆に2^nなら前半だ。次にt/2の配列で考えて、それが2^(n-1)+1と表せたら(すなわち、t/2が奇数だったら、さらに言えばtの1ビット目が1だったら=>t&2が1だったら)、その配列の後半に来るはずだ。ここまでで分かったのは、F(t)が入る場所は、(t&1)*2^(n-1) + (t&2)*2^(n-2)+ 「何か数」という式で表すことができる。これをずっと続けて
(t&1)*2^(n-1) + (t&2)*2^(n-2) + (t&4)*2^(n-3)………+(t&(2^(n-1)))*2(n-n)
という式で与えられる。


ビット列を扱うことに長けている方ならティーンとしたことだろう。
そう、tをnbitで表現した場合、tのビット列を反転した値になっているのだ。これを「ビットの反転」と呼ぶ。



なお、ビットの反転アルゴリズムの解説までは骨が折れるので、FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法を参照していただきたい。ここでは、「リスト1.2.5-3. より効率的なビット反転並べ替え」をまるっとそのまま引用した。



このソースコードは次のようになる。
なお、Javaは配列からサブ配列が取り出せない糞仕様なので、配列の先頭番地(arrayIndex)を渡すことで解決した。


ソース2

public static void fft(double[] real,double[] image){
  int n = real.length;
  double theta = -2*Math.PI/n;
  _fft(n,theta,real,image,0);

  //並べ替え
  int j,i,N=n,nh1,nh;
  for (j = 0,i=0,nh = N>>1,nh1 =nh+1; j < nh; j += 2) {
    double tmpr,tmpi;
      if (j < i) {
          tmpr = real[j];
          real[j] = real[i];
          real[i] = tmpr;
          tmpr = real[j + nh1];
          real[j + nh1] = real[i + nh1];
          real[i + nh1] = tmpr;
          tmpi = image[j];
          image[j] = image[i];
          image[i] = tmpi;
          tmpi = image[j + nh1];
          image[j + nh1] = image[i + nh1];
          image[i + nh1] = tmpi;
      }
      tmpr = real[j + nh];
      real[j + nh] = real[i + 1];
      real[i + 1] = tmpr;

      tmpi = image[j + nh];
      image[j + nh] = image[i + 1];
      image[i + 1] = tmpi;
      for (int k = nh >> 1; k > (i ^= k); k >>= 1);
  }
}

private static void _fft(
    int n,double theta,
    double[] real,double[] image,int arrayIndex){
  if(n<=1)return;
  int nh = n/2;
  for(int m=0;m<nh;m++){
    double r1 = real[m+arrayIndex],r2 = real[m+nh+arrayIndex];
    double i1 = image[m+arrayIndex],i2 = image[m+nh+arrayIndex];
    real[m+arrayIndex] = r1+r2;
    image[m+arrayIndex] = i1+i2;

    double Rm = r1 - r2;
    double Im =  i1 - i2;
    double cos = Math.cos(m*theta);
    double sin = Math.sin(m*theta);
    real[m+nh+arrayIndex] = Rm*cos-Im*sin;
    image[m+nh+arrayIndex] = Rm*sin+Im*cos;
  }
  _fft(nh,2*theta,real,image,arrayIndex);
  _fft(nh,2*theta,real,image,arrayIndex+nh);
}

先と同じ条件でのこれの実行時間は1.8秒程度。ちょっと速くなったね。何よりもメモリの使用がぐっと改善されました。
次回、より速く動作させるための改善を行っていきます。