复数(Complex)类及FFT的C++实现

===○专业造轮子○===

因为一些原因需要额外编写傅里叶变换(Fourier Transform)的实现代码，

而傅里叶变换需要复数的支持，因此额外编写了一个复数类。

首先是复数(Complex)类的设计，很简单，数据域只有实部和虚部，为了方便外部操作，

将数据设置为公开访问，大致如下

classComplex

{

public:

    Complex();

    Complex(double re,double im);

   
// operator +,-,*,/,= etc.

   
double real;  // 实部

   
double imag;  // 虚部

};

为了简便操作，Complex类重载了基本运算符。

接下来是傅里叶变换的实现。

为了提高效率，在实际应用中(离散形式的傅里叶变换,DFT)多采用快速算法即FFT

全称为Fast Fourier Transform

此处引用到复数类Complex,稍后详述。

FFT的具体实现如下

函数声明FFT()

#ifndef FOURIER_H
#define FOURIER_H

#include "Complex.h"

#ifndef CONST_PI
#define CONST_PI
const double PI = 3.14159265358979;
const double PI_X2 = 2 * PI;
#endif

/* r=log2(N) */
extern void FFT(Complex *TD, Complex*FD, int r);

#endif

函数实现FFT()

#include "Fourier.h"
#include <stdlib.h>
#include <math.h>

/* r=log2(N) */
void FFT(Complex *TD, Complex*FD, int r)
{
const int count = 1 << r;
int csz = sizeof(Complex)*count;
Complex* W = (Complex*)malloc(csz / 2);
Complex* X1 = (Complex*)malloc(csz);
Complex* X2 = (Complex*)malloc(csz);
Complex* X = NULL;

int  i, j, k;
int  dist, p;
double f = PI_X2 / count;
double a = 0;
for (i = 0; i < count / 2; ++i)
{
W[i] = Complex(cos(a), -sin(a));
a += f;
}

for (i = 0; i < count; ++i)
{
X1[i] = TD[i];
}

for (k = 0; k < r; ++k)
{
for (j = 0; j < (1 << k); ++j)
{
dist = 1 << (r - k);
for (i = 0; i < dist / 2; ++i)
{
p = j * dist;
X2[i + p] = X1[i + p] + X1[i + p + dist / 2];
X2[i + p + dist / 2] = (X1[i + p] - X1[i + p + dist / 2])* W[i * (1 << k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}

for (j = 0; j < count; ++j)
{
p = 0;
for (i = 0; i < r; ++i)
{
if (j&(1 << i))
{
p += 1 << (r - i - 1);
}
}
FD[j] = X1[p];
}

free(W);
free(X1);
free(X2);
}

复数类Complex的声明文件Complex.h

#ifndef COMPLEX_H
#define COMPLEX_H

#ifndef BOOLEAN_VAL
#define BOOLEAN_VAL
typedef int BOOL;
#define TRUE 1
#define FALSE 0
#endif

class Complex
{
public:
Complex();
Complex(double re, double im);

Complex operator=(double v);
Complex operator+(double v);
Complex operator-(double v);
Complex operator*(double v);
Complex operator/(double v);
Complex operator+=(double v);
Complex operator-=(double v);
Complex operator*=(double v);
Complex operator/=(double v);

Complex operator=(Complex c);
Complex operator+(Complex c);
Complex operator-(Complex c);
Complex operator*(Complex c);
Complex operator/(Complex c);
Complex operator+=(Complex c);
Complex operator-=(Complex c);
Complex operator*=(Complex c);
Complex operator/=(Complex c);

BOOL operator==(Complex c);
BOOL operator!=(Complex c);

double real;
double imag;
};

#endif

复数类Complex的实现文件Complex.cpp

#include "Complex.h"

Complex::Complex()
{
real = 0;
imag = 0;
}

Complex::Complex(double re, double im)
{
real = re;
imag = im;
}

Complex Complex::operator+(double v)
{
return Complex(real + v, imag);
}

Complex Complex::operator-(double v)
{
return Complex(real - v, imag);
}

Complex Complex::operator*(double v)
{
return Complex(real*v, imag*v);
}

Complex Complex::operator/(double v)
{
return Complex(real / v, imag / v);
}

Complex Complex::operator=(double v)
{
real = v;
imag = 0;
return *this;
}

Complex Complex::operator+=(double v)
{
real += v;
return *this;
}

Complex Complex::operator-=(double v)
{
real -= v;
return *this;
}

Complex Complex::operator*=(double v)
{
real *= v;
imag *= v;
return *this;
}

Complex Complex::operator/=(double v)
{
real /= 2;
imag /= 2;
return *this;
}

Complex Complex::operator+(Complex c)
{
return Complex(real + c.real, imag + c.imag);
}

Complex Complex::operator-(Complex c)
{
return Complex(real - c.real, imag - c.imag);
}

Complex Complex::operator*(Complex c)
{
double re = real*c.real - imag*c.imag;
double im = real*c.imag + imag*c.real;
return Complex(re, im);
}

Complex Complex::operator/(Complex c)
{
double x = c.real;
double y = c.imag;
double f = x*x + y*y;
double re = (real*x + imag*y) / f;
double im = (imag*x - real*y) / f;
return Complex(re, im);
}

Complex Complex::operator=(Complex c)
{
real = c.real;
imag = c.imag;
return *this;
}

Complex Complex::operator+=(Complex c)
{
real += c.real;
imag += c.imag;
return *this;
}

Complex Complex::operator-=(Complex c)
{
real -= c.real;
imag -= c.imag;
return *this;
}

Complex Complex::operator*=(Complex c)
{
double re = real*c.real - imag*c.imag;
double im = real*c.imag + imag*c.real;
real = re;
imag = im;
return *this;
}

Complex Complex::operator/=(Complex c)
{
double x = c.real;
double y = c.imag;
double f = x*x + y*y;
double re = (real*x + imag*y) / f;
double im = (imag*x - real*y) / f;
real = re;
imag = im;
return *this;
}

BOOL Complex::operator==(Complex c)
{
return ((real == c.real) && (imag == c.imag));
}

BOOL Complex::operator!=(Complex c)
{
return ((real != c.real) || (imag != c.imag));
}

在VS2013中编写测试程序，以16384个复数组成的数组为样本进行测试，从启动到输出结果远不足1秒，证明FFT效率较高。

以同样的数据在Matlab中进行测试，发现结果一致，由此证明以上算法代码基本无误(个人编辑代码可能会出现小的疏忽)。



本文原创，博文原始地址

http://blog.csdn.net/fengyhack/article/details/42780785