# | Likes | Tech tags | Title | Creator | Created date |
---|---|---|---|---|---|
1 | 0 |
2022-11-25 22:19
|
A commented and somewhat optimized iterative implementation of the Cooley-Tukey fast Fourier transform algorithm. Runs in \(O(n \log n)\) time and \(O(n)\) space (where \(n\) is the size of the input array).
Important note: The size of the input array must be a power of two
#define _USE_MATH_DEFINES
#include <complex>
#include <vector>
#include <cmath>
// Assumes the size of x is a power of 2
std::vector<std::complex<double> > fast_fourier_transform(std::vector<std::complex<double> > x, bool inverse = false) {
std::vector<std::complex<double> > w(x.size(), 0.0); // w[k] = e^{-i2\pi k/N}
// Precalculate w[k] for faster FFT computation, do it in a 'binary-search' way to provide decent numeric accuracy
w[0] = 1.0;
for(int pow_2 = 1; pow_2 < (int)x.size(); pow_2 *= 2) {
w[pow_2] = std::polar(1.0, 2*M_PI * pow_2/x.size() * (inverse ? 1 : -1) );
}
for(int i=3, last=2; i < (int)x.size(); i++) {
// This way of computing w[k] guarantees that each w[k] is computed in at most log_2 N multiplications
if(w[i] == 0.0) {
w[i] = w[last] * w[i-last];
} else {
last = i;
}
}
for(int block_size = x.size(); block_size > 1; block_size /= 2) {
// Do the rearrangement for 'recursive' call block by block for each level'
std::vector<std::complex<double> > new_x(x.size());
for(int start = 0; start < (int)x.size(); start += block_size) {
for(int i=0; i<block_size; i++) {
new_x[start + block_size/2 * (i%2) + i/2] = x[start + i];
}
}
x = new_x;
}
for(int block_size = 2; block_size <= (int)x.size(); block_size *= 2) {
// Now compute the FFT 'recursively' level by level bottom-up
std::vector<std::complex<double> > new_x(x.size());
int w_base_i = x.size() / block_size; // w[w_base_i] is the e^{-i2\pi / N} value for the N of this level
for(int start = 0; start < (int)x.size(); start += block_size) {
for(int i=0; i < block_size/2; i++) {
new_x[start+i] = x[start+i] + w[w_base_i*i] * x[start + block_size/2 + i];
new_x[start+block_size/2+i] = x[start+i] - w[w_base_i*i] * x[start + block_size/2 + i];
}
}
x = new_x;
}
return x;
}
int main() {
std::vector<std::complex<double> > result = fast_fourier_transform({1, 4, 3, 6, 2, 3, 7, 4});
return 0;
}
functions | ||
std::polar |
en.cppreference.com | cplusplus.com |
std::vector::operator[] |
en.cppreference.com | cplusplus.com |
std::vector::size |
en.cppreference.com | cplusplus.com |
std::vector::vector |
en.cppreference.com | cplusplus.com |
classes | ||
std::complex |
en.cppreference.com | cplusplus.com |
std::vector |
en.cppreference.com | cplusplus.com |
Implement the Cooley–Tukey algorithm for computing fast Fourier transform. Use the implemented algorithm to compute the FFT for some set of inputs.