高速フーリエ変換のソースコードです。 「高速フーリエ変換(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; } |