Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Need help understanding this line in an FFT algorithm

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 twiddle factor. 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 another twiddle factor, and I can see that equivalence, 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?

like image 718
Verpous Avatar asked Oct 05 '20 18:10

Verpous


People also ask

What is the FFT algorithm?

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:

How many lines of code does FFT have?

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!

What is the fast Fourier transform (FFT)?

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.

How do you calculate an FFT?

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.


1 Answers

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.

like image 138
harold Avatar answered Oct 07 '22 02:10

harold