FFT (高速フーリエ変換) で信号の高周波数の部分をカットする

信号処理は専門じゃない(どころかそういう授業を受けたりするような環境ですら無かった)ので全然わからないが,ひょんなことからちょっと信号を処理することになった.

FFT は普段からたまに使うのだけれど,いつもただ畳込みを計算するばかりで,周波数がどうのみたいなことをやったことがなかった.以下がプログラミングコンテストでいつも使っている FFT のコード.(a の長さは 2 の累乗にする.)

#include <complex>
typedef complex<double> c_t;

void fft(vector<c_t> &a, int dir = 1) {
  int n = a.size(), i = 0;
  rep (j, n) {
    if (j < i) swap(a[i], a[j]);
    for (int k = n >> 1; k > (i ^= k); k >>= 1);
  }

  for (int s = 0; 1 << s < n; s++) {
    int m = 1 << s;
    c_t wm = polar(1.0, M_PI / m * dir), w(1.0);
    rep (j, m) {
      for (int k = j; k < n; k += m * 2) {
        c_t s = a[k], t = w * a[k + m];
        a[k] = s + t;
        a[k + m] = s - t;
      }
      w *= wm;
    }
  }
}

これにぶち込んだあと,配列のある部分を 0 にして,またぶち込む(dir = -1 にして呼び出し,さらに全てを a.length() で割り,実数部を取る)と,高周波数領域がカットできる.

これでどこを 0 にすればよいかというと,配列の真ん中の部分である.これはちょっと意外な感じがする(最初の部分とか後半の部分というのがよくある感じがする)が,まぁよく考えたらそうである.