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 にすればよいかというと，配列の真ん中の部分である．これはちょっと意外な感じがする（最初の部分とか後半の部分というのがよくある感じがする）が，まぁよく考えたらそうである．