In my program I have a function that performs the fast Fourier transform. I know there are very good implementations freely available, but this is a learning thing so I don't want to use those. I ended up finding this comment with the following implementation (it originated from the Italian entry for the FFT):
void transform(complex<double>* f, int N) //
{
ordina(f, N); //first: reverse order
complex<double> *W;
W = (complex<double> *)malloc(N / 2 * sizeof(complex<double>));
W[1] = polar(1., -2. * M_PI / N);
W[0] = 1;
for(int i = 2; i < N / 2; i++)
W[i] = pow(W[1], i);
int n = 1;
int a = N / 2;
for(int j = 0; j < log2(N); j++) {
for(int k = 0; k < N; k++) {
if(!(k & n)) {
complex<double> temp = f[k];
complex<double> Temp = W[(k * a) % (n * a)] * f[k + n];
f[k] = temp + Temp;
f[k + n] = temp - Temp;
}
}
n *= 2;
a = a / 2;
}
free(W);
}
I've made a lot of changes by now but this was my starting point. One of the changes I made was to not cache the twiddle factors, because I decided to see if it's needed first. Now I've decided I do want to cache them. The way this implementation seems to do it is it has this array W
of length N/2
, where every index k
has the value . What I don't understand is this expression:
W[(k * a) % (n * a)]
Note that n * a
is always equal to N/2
. I get that this is supposed to be equal to , and I can see that , which this relies on. I also get that modulo can be used here because the twiddle factors are cyclic. But there's one thing I don't get: this is a length-N DFT, and yet only N/2
twiddle factors are ever calculated. Shouldn't the array be of length N, and the modulo should be by N?
The FFT is a fast, O[NlogN] algorithm to compute the Discrete Fourier Transform (DFT), which naively is an O[N2] computation. The DFT, like the more familiar continuous version of the Fourier transform, has a forward and inverse form which are defined as follows:
Use that. The FFT routines here have less than a hundred lines of code. The library implements forward and inverse fast Fourier transform (FFT) algorithms using both decimation in time (DIT) and decimation in frequency (DIF). @DaBler That's exactly what I was searching for! thank you!
One such method was developed in 1965 by James W. Cooley and John W. Tukey 1 Their work led to the development of a program known as the fast Fourier transform. The fast Fourier transform (FFT) is a computationally efficient method of generating a Fourier transform.
To calculate an FFT (Fast Fourier Transform), just listen. The human ear automatically and involuntarily performs a calculation that takes the intellect years of mathematical education to accomplish.
But there's one thing I don't get: this is a length-N DFT, and yet only N/2 twiddle factors are ever calculated. Shouldn't the array be of length N, and the modulo should be by N?
The twiddle factors are equally spaced points on the unit circle, and there is an even number of points because N is a power-of-two. After going around half of the circle (starting at 1, going counter clockwise above the X-axis), the second half is a repeat of the first half but this time it's below the X-axis (the points can be reflected through the origin). That is why Temp
is subtracted the second time. That subtraction is the negation of the twiddle factor.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With