One-dimensional fast Fourier transform code

Source: Internet
Author: User
Tags setf

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

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.