Quick C++ implementation of Fast Fourier Transform

Content

1. Background and motivation

2. Literation

3. Butterfly Computing based on literation

Released Version

---

1. Background and motivation

As one of the greatest and most profound mathematical discoveries of all time, the Fourier Transform is widely and deeply used in mathematics and physics, especially in the field of communication engineering.

Discrete Fourier transform, DFT, is the most basic method of signal analysis, through which the signal is transformed from the time domain to the frequency domain, and then the spectrum structure and change rule of the signal is studied.

However, for a complex input signal, for example, a complex sequence of N terms, X[n] , it seems that it needs about N*N times of calculation to obtain the output signal. If N is a large number, the process of DFT is obviously a really huge project for user’s computer. In order to preserve memory and transform input signal more efficiently, the code ultimately does Fast Fourier Transform (FFT), which is a fast algorithm of DFT.

FFT can significantly reduce the computation requirements. Generally speaking, FFT uses symmetry and periodicity of wn to halving the amount of terms of a sequence, thereby only (N/2)2 times of computation are needed for each part. As number of halving increases, then computation requirements rapidly decrease. The more halving, the greater the savings, which is the advantage of FFT.

---

2. Literation

Literation of input array, whose elements are sampling points, is an appropriate and efficient method. The main idea is that:

‘Regard each TTT as a separate polynomial, then combine them two by two to form a polynomial (each polynomial has 2 terms), and continue to combine (the second polynomial has 4 terms each) until only one polynomial (with N terms) is left.’

Such a code implementation of this algorithm is simple and memory-friendly as this method can complete the calculation of multiple polynomials with the least number of computations and then save the memory as possible.

3. Butterfly Computing based on literation

After we understand how to reduce the number of computations through literation, now we need to achieve the calculation of coefficients.

The key is to use the symmetry and periodicity of DFT coefficients.

The advantage of butterfly computing based on literation is that traditional implementation of Butterfly Processing is based on that several for loops are nested inside each other. However, my original method to achieve this algorithm is function’s literation. Only one for loop is used in FFT function. Although there's not much reduction in computation, the code is more compact and readable.

A whole process is shown below:

1.Classification

Divide x(n) into even and odd two groups according to n:

When n is even, n = 2r.

When n is odd, n=2r+1.

Then we get

Code:

2.Literation

We need do classification over and over again until it cannot continue to divide.

Code:

3.Calculation

Through proper simplification by using the symmetry and periodicity of DFT coefficients, we get:

 

Then a Butterfly Processing Element is generated.

The instructions for this butterfly processing element are as follows:

1. One butterfly operation requires one multiplication and two multiplications.
2. The two on the left are inputs.
3. The two on the right are outputs.
4. A small circle in the middle represents the addition and subtraction operations (the path above is the addition output, the path below is the subtraction output).

For example, when N = 8, the butterfly process will be :

The above three procesures are complete FFT calculation process.

---

Released Version

#include <complex>
#include <iostream>
#include <valarray>

using namespace std;
 
const double PI = 3.141592653589793238460;
 
typedef std::complex<double> Complex;
typedef std::valarray<Complex> ComplexArray;
 
// FFT
void fft(ComplexArray &a)
{
	const size_t N = a.size();
	if (N <= 1) return;

	// classification
	ComplexArray even = a[std::slice(0, N/2, 2)];	 //+ 
	ComplexArray odd  = a[std::slice(1, N/2, 2)];	 //-
 
	// literation
	fft(even);
	fft(odd);
 
	// calculaiton
	for (size_t k = 0; k < N/2; k++)
	{
		Complex rotating = std::polar(1.0, -2 * PI * k/N) * odd[k];
		a[k] = even[k] + rotating;
		a[k + N / 2] = even[k] - rotating;
	}
}


// inverse FFT
void ifft(ComplexArray& a)
{
	// conjugate complex numbers
	a = a.apply(std::conj);
 
	// quick literation
	fft(a);
 
	// conjugate complex numbers
	a = a.apply(std::conj);
 
	// 1/N
	a /= a.size();
}
 
int main()
{
	const int N=8;
    Complex input[N]; 
    cout<<"Please input:\n";

    for(int i=0;i<N;i++){

       cout<<"input["<<i<<"]=";
       cin>>input[i];
    }

	ComplexArray sample_points(input, N);
 
	fft(sample_points);
 
	std::cout << "fft" << std::endl;
	for (int i = 0; i < N; i++){

		std::cout << sample_points[i] << std::endl;
	}
 
	ifft(sample_points);
 
	std::cout << std::endl << "ifft" << std::endl;
	for (int i = 0; i < N; i++){

		std::cout << sample_points[i] << std::endl;
	}
	
	return 0;
}