高速フーリエ変換のソースコードです。 「高速フーリエ変換(FFT)をおじさんもC++で作ってみたよ」を参考にしました。

void fft(vector< complex<double> > *pin, vector<complex< double> > *pout){
  int n_level;
  int size = pin->size();

  int i = 0;
  for(i = 0; i < 64; i++){
    if(size >> i == 1){
      break;
    }
  }

  n_level = i;

  vector<int> id;

  id.push_back(0);
  id.push_back(1);

  for(int i = 0; i < n_level - 1; i++){
    int id_size = id.size();
    for(int j = 0; j < id.size(); j++){
      id[j] = id[j] * 2;
    }

    id.insert(id.end(), id.begin(), id.end());
    int k = id_size;
    for(k = id_size; k < id.size(); k++){
      id[k] = id[k] + 1;
    }
  }

  pout->resize(size);
  for(int i = 0; i < size; i++){
    (*pout)[i] = (*pin)[id[i]];
  }

  unsigned int po2 = 1;
  for(int i_level = 1; i_level <= n_level; i_level++){
    po2 <<= 1;
    int po2m = po2 >> 1;

    complex<double> w = exp(std::complex<double>(0.0, 2 * M_PI/(double)po2));
    complex<double> ws = complex<double>(1, 0);

    for(int k = 0; k < po2m; k++){
      for(int j = 0; j < size; j += po2){
        complex<double> *pa = &(*pout)[j + k];
        complex<double> *pb = &(*pout)[j + k + po2m];
        complex<double> wfb = ws * *pb;
        *pb = *pa - wfb;
        *pa += wfb;
      }
      ws *= w;
    }
  }

  return;
}

  添付編集
Last-modified: 2014-08-01 (金) 20:49:15 (3527d)