Fast Fourier Transform C ++ Recursive Algorithm Implementation

Source: Internet
Author: User

Fast Fourier Transform C ++ Recursive Algorithm Implementation

 

Some algorithms on the Internet are incorrectly tested and run, although the Code is non-recursive. To verify the accuracy of the fast Fourier transform, I provide a self-designed recursive algorithm.

Code of the "base 2" Fast Fourier Transform Algorithm Based on Time Domain extraction:

Fouier. h file:

# Pragma once
# Include "complex. H"
Class Fouier
{
Complex * data;
Void FFT (INT start, int step, int Len );
Complex W (int K, int N); // e ^ (-I * 2 * pI * k/N)
Public:
Fouier (void );
~ Fouier (void );

Void FFT ();
};

Fouier. c file:

# Include "Fouier. H"
# Include <iostream>
Using namespace STD;
# Include <cmath>
# Include <ctime>
# Define datalen 32
# Define keyValue 10000 // generate a random floating point value to ensure that the numerator and denominator are within this value
# Define PI 3.14159265358979323846
Fouier: Fouier (void)
{
Data = new complex [datalen];
Srand (unsigned int (time (0 )));
Cout <"source data:" <Endl;
For (INT I = 0; I <datalen; I ++)
{
Data [I] = (RAND () % (keyValue)/(double) (RAND () % (keyValue) + 1 );
If (I % 5 = 0 & I! = 0)
Cout <Endl;
Cout <data [I] <"";
}
Cout <Endl;
}

Fouier ::~ Fouier (void)
{
Delete [] data;
}
Complex Fouier: w (int K, int N) // Euler's Formula
{
Double alpha =-2 * pI * k/N;
Return complex (COS (alpha), sin (alpha ));
}
Void Fouier: FFT (INT start, int step, int Len)
{
If (LEN = 1) // an element
{
// One element does not need to be transformed.
Return;
}
FFT (START, step * 2, Len/2); // x1 (k)
FFT (start + step, step * 2, Len/2); // X2 (k)
Complex x1, x2;
For (INT I = 0; I <Len/2; I ++)
{
X1 = data [start + step * I * 2];
X2 = data [start + step * (I * 2 + 1)];
// Calculate X (k): K = 0 ~ N/2-1
Data [start + step * I] = X1 + W (I, Len) * X2;
// Calculate X (k): K = n/2 ~ N-1
Data [start + step * (I + Len/2)] = X1-W (I, Len) * X2;
}
}
Void Fouier: FFT ()
{
FFT (0, 1, datalen );
Cout <"transformed data:" <Endl;
For (INT I = 0; I <datalen; I ++)
{
If (I % 5 = 0 & I! = 0)
Cout <Endl;
Cout <data [I] <"";
}
}

Complex. h file:

# Pragma once
# Include <iostream>
Using namespace STD;
Class complex // A + B * I
{
Double A; // real part
Double B; // virtual part
Public:
Complex (double A = 0, double B = 0 );
// + Operation
Friend complex operator + (complex & X, complex & Y );
Friend complex operator + (Double X, complex & Y );
Friend complex operator + (complex & X, Double Y );
//-Operation
Friend complex operator-(complex & X, complex & Y );
Friend complex operator-(Double X, complex & Y );
Friend complex operator-(complex & X, Double Y );
// * Operation
Friend complex operator * (complex & X, complex & Y );
Friend complex operator * (Double X, complex & Y );
Friend complex operator * (complex & X, Double Y );
// = Operation
Complex operator = (complex & X );
Complex operator = (Double X );
// <Operation
Friend ostream & operator <(ostream & out, complex & C );

~ Complex (void );
};

Complex. c file:

# Include "complex. H"

Complex: complex (double A, double B) // The default virtual part is 0.
{
This-> A =;
This-> B = B;
}

Complex ::~ Complex (void)
{
}

Complex operator + (complex & X, complex & Y)
{
Return complex (X. A + Y. A, X. B + Y. B );
}
Complex operator + (Double X, complex & Y)
{
Return complex (X + Y. A, Y. B );
}
Complex operator + (complex & X, Double Y)
{
Return complex (X. A + Y, X. B );
}

Complex operator-(complex & X, complex & Y)
{
Return complex (X. a-y.a, X. b-y. B );
}
Complex operator-(Double X, complex & Y)
{
Return complex (x-y.a,-Y. B );
}
Complex operator-(complex & X, Double Y)
{
Return complex (X. A-y, X. B );
}

Complex operator * (complex & X, complex & Y)
{
Return complex (X. A * Y. a-x. B * Y. B, X. A * Y. B + X. B * Y. );
}
Complex operator * (Double X, complex & Y)
{
Return complex (x * Y. A, x * Y. B );
}
Complex operator * (complex & X, Double Y)
{
Return complex (X. A * y, X. B * y );
}

Complex complex: Operator = (complex & X)
{
A = X.;
B = x. B;
Return * this;
}
Complex complex: Operator = (Double X)
{
A = X;
B = 0;
Return * this;
}
Ostream & operator <(ostream & out, complex & C)
{
If (C.! = 0 | C. A = 0 & C. B = 0)
Out <c.;
If (C. B! = 0)
{
If (C. B> 0)
Out <"+ ";
If (C. B! = 1)
Out <C. B;
Out <"I ";
}
Return out;
}

Main. c file:

# Include <iostream>
Using namespace STD;
# Include "Fouier. H"

Int main ()
{
Fouier F;
F. FFT ();
Return 0;
}

If you have any mistakes, please correct them!

References: http://zhoufazhe2008.blog.163.com/blog/static/63326869200971010421361/

Wikipedia-http://zh.wikipedia.org/wiki/%E5%BF% AB %E9%80%9F%E5%82%85%E7% AB %8B%E5%8F%B6%E5%8F%98%E6%8D%A2

Related Article

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.