The previous essay briefly wrote an algorithm for the reordering of arrays in FFT. Now share the full FFT code (with more detailed comments).
/*2015 November 10 at Hebei University of Technology * *
#include <complex>
#include <iostream.h>
#include <math.h>
#include <stdlib.h>
const int n=8; The length of the array
const double pi=3.141592653589793; Pi
const double zero= 1.0E-14; When the absolute value is less than this number, it is 0.
typedef std::complex<double> Complex;
Function Prototype Declaration
void Reverse (complex src[],int N);
void FFT (complex src[],int N);
void Display (complex f[],int N);
/* Main function */
int main ()
{
Complex src[8]; Original array
for (int i=0;i<8;i++)
{
Src[i]=complex (double (i), 0.0);
}
cout<< "original array:" <<endl;
Display (Src,n);
One-dimensional fast Fourier transform
FFT (Src,n);
cout<< "transformed array is:" <<endl;
Display (Src,n);
return 0;
}
/***********************************************
One-dimensional fast discrete Fourier transform *************
Input: Array length n *************
Imaginary part imag *************
Real part real *************
************************************************/
void FFT (complex src[],int N)
{
Original array in descending order of bit code
Reverse (Src,n);
cout<< "Inverted array is:" <<endl;
Display (Src,n);
Calculated M,m=log (2,n) indicates a total M-level butterfly operation
for (int i=n,m=1; (I=I/2)!=1; m++);
for (int l=1; l<=m; l++)
{
int Arraynum=n/pow (2,L); The number of small arrays that are divided
int Elementnum=pow (2,L); The number of elements per decimal group
int UNITEACHARR=ELEMENTNUM/2; The number of butterfly units included in each decimal group
int Offset=pow (2,L-1); The distance of each butterfly operation node
cout<< "Offset:" <<offset<<endl;
for (int i=0;i<n;i=i+elementnum)
{
for (int j=i;j<=i+uniteacharr-1;j++)
{
Double r= (j-i) *pow (2,m-l); R in W (R,n)
Complex W=complex (cos (2*pi*r/n), sin (2*pi*r/n));
Complex Temp=src[j+offset]*complex (cos (2*pi*r/double (n)),-sin (2*pi*r/double (n)));
Src[j+offset]=src[j]-temp;
Src[j]=src[j]+temp;
}
}
}
}
/***********************************************
function function: Implement original array reordering *******
The main implementation of the position code in reverse order * * * *
Function Input: Original array * * * *
************************************************/
void Reverse (complex src[],int N)
{
/*
P1 and P2 represent the index of the two elements to be exchanged
Where p1 from the second element to the second element of the penultimate
Scanning, p2 is determined by the algorithm and P1
The index of another element of the interchange
*/
int p1;
int p2;
int middle; Represents a center point
int offset; Represents an offset
/* Center point and offset are important parameters to calculate P2 */
P2=N/2; The index of the element exchanged with P1=1 is N/2
For (p1=1;p1<=n-2;/* first element and last element without swapping */p1++)
{
/* Swap when P1<P2, avoid duplicate exchange */
if (P1<P2)
{
/* Swap */
Complex TEMP=SRC[P1];
SRC[P1]=SRC[P2];
Src[p2]=temp;
}
/* Start calculating next p2*/below
MIDDLE=N/2; The initial center point
OFFSET=P2; The offset of the preliminary examination
To update when the center point is less than or equal to the offset
while (Middle<=offset)
{
Offset=offset-middle;
MIDDLE=MIDDLE/2;
}
The P2=MIDDLE+OFFSET;//P2 is equal to the new center point and the offset and
}
}
/* Show Data */
void Display (complex f[],int N)
{
for (int i=0;i<n;i++)
{
Cout.width (9);
COUT.SETF (ios::fixed);
Cout.precision (4);
Cout<<f[i].real ();
Cout.flags (Ios::showpos);
Cout.width (9);
COUT.SETF (ios::fixed);
Cout.precision (4);
Cout<<f[i].imag () << ' i ' << ' t ';
COUT.UNSETF (Ios::showpos);
if ((i+1)%3==0) cout<<endl;
}
cout<<endl;
}
One-dimensional fast Fourier transform code