C + + implementation of fast Fourier transform (FFT) and MATLAB experiment

By Jia Jia's "complex transformation function and integral transform" for two days, finally understand how the Fourier transform is the same thing. But to achieve fast Fourier transform but do not need to understand so many things, see the "Introduction to the algorithm" in the 30th chapter "polynomial and Fast Fourier transform" on it. However, the introduction of algorithmic introduction and the standard of a little different, that is, the rotation factor is just the reverse, but still equivalent.

Standard discrete Fourier transform forms such as:

Yk=σj=0n-1 ajωn-kj = A (ωn-k).

(Ωnk is the k n-th square root of complex number 1 and defines polynomial A (x) =σj=0n-1 AJXJ)

The discrete Fourier inverse transformation idft (inverse DFT) forms such as: aj= (σk=0n-1 ykωnkj)/N.

The following different color contents are referenced and amended:

Fast Fourier transform (FAST FOURIER TRANSFORM,FFT) is a fast algorithm for discrete Fourier transforms (discrete Fourier transform,dft), which is based on discrete Fourier transforms of singular, even, virtual, Real property, the algorithm of discrete Fourier transform is improved. It has no new discoveries about the theory of Fourier transform, but it is a big step to apply the discrete Fourier transform to the computer system or the digital systems. The
sets the complex sequence of Xn to N items, by DFT, the computation of any Xi requires N times plural multiplication and N-1 plural addition, while a complex multiplication is equal to four real multiplication and two real addition, one complex addition equals two real addition, even if a plural multiplication and a plural addition are defined as a " Operation "(four times real multiplication and four real addition), then it is necessary to find out the Xi of n-order complex sequence, that is, the N-point DFT transformation requires N2. When N = 1024 points or more, requires N2 = 1,048,576 operations, in FFT, using ωn periodicity and symmetry, an n-item sequence (set n is even), divided into two n/2 items of the sequence, each N/2-point DFT Transform needs (n/ 2) 2 times operation, and then using N-time operation to combine the DfT transform of two n/2 points into a DfT transformation of n points. After this transformation, the total number of operations becomes n + 2 * (N/2) 2 = n + N2/2. Continuing with the example above, when N =1024, the total number of operations becomes 525,312 times, saving about  50% operations. And if we keep this "split" thinking going, until divided into 221 groups of DFT operation Unit, then the N-point DFT transformation requires only n * log2n times, n = 1024 points, the operation is only 10,240 times, is the previous direct algorithm 1%, the more points, the operation The greater the amount of savings, this is the superiority of the FFT.

The implementation of FFT can be top-down, recursive, but for hardware to achieve high cost, software implementation is not efficient, instead of iterative better, from the bottom to solve the problem. The iterative version of the sense and merge sort is similar, but first the "bit reversal permutation" method is used to put Xi in the right place, set I and J to each other s = log2n bit binary palindrome number, assuming s = 3, I = (2 = 6, J = (011) 2 = 3, if I≠j, then Xi and Xj should swap positions. (The generation of this palindrome number, is very interesting and is very basic operation, want to have a beginner in C + + when there is such a problem.) When the "bit reversal permutation" is complete, consider each XI as a separate polynomial, and then merge them into a polynomial (2 for each polynomial) in two or two, and the merging is actually "Butterfly" (butterfly Operation, refer to the introduction to algorithms ^_^), Continue merging (4 for each polynomial for the second time) until only one polynomial (with N) is left, so that the merged layer is LOG2N, each layer has n operations, so there is a total of n * log2n operations. The iterative process, as shown in the following illustration, is bottom-up.

Write your own C + + source program, as follows:

/**/////////////////////////////////////////////////
//
Rapid Fourier transform fast Fourier Transform
by rappizit@yahoo.com.cn
2007-07-20
Version 2.0
Improved the algorithm introduction of algorithm, rotation factor to take ωn-kj (ωnkj conjugate plural)
And only N/2 times are calculated, and the need to compute (n * LG N)/2 times is not improved.
//
/**/////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const int N = 1024;
const float PI = 3.1416;

inline void swap (float &a, float &b)
... {
float T;
t = A;
A = b;
b = t;
}

void Bitrp (float xreal [], float ximag [], int n)
... {
Bit reversal permutation bit-reversal permutation
int I, J, A, B, p;

for (i = 1, p = 0; i < n; I *= 2)
... {
p + +;
}
for (i = 0; i < n; i + +)
... {
A = i;
b = 0;
for (j = 0; J < P; j + +)
... {
b = (b << 1) + (A & 1); b = b * 2 + a% 2;
a >>= 1; A = A/2;
}
if (b > i)
... {
Swap (Xreal [i], xreal [b]);
Swap (Ximag [i], Ximag [b]);
}
}
}

void FFT (float xreal [], float ximag [], int n)
... {
Fast Fourier transform, the complex X transformation is still preserved in X, Xreal, Ximag is the real part of X and imaginary part
float wreal [N/2], Wimag [N/2], Treal, Timag, Ureal, u