?好象不斷有人問關(guān)于FFT的問題，我也發(fā)一個(gè)以前寫的FFT程序吧。先發(fā)一個(gè)復(fù)數(shù)運(yùn)算的，過幾天再發(fā)一個(gè)定點(diǎn)的。我自己覺得這個(gè)程序的可讀性還是比較好的，如果你能理解蝶形運(yùn)算，你就應(yīng)該能理解這個(gè)程序。我的程序沒有用通用的STL復(fù)數(shù)類，而用了我以前自己寫的一個(gè)復(fù)數(shù)運(yùn)算類，只是做了一些STL化，大學(xué)里的老師應(yīng)該不會(huì)罵我了吧:)。我的這個(gè)復(fù)數(shù)類也有其優(yōu)點(diǎn)，可以支持直角坐標(biāo)和極坐標(biāo)的交互運(yùn)算。比如說，假設(shè)第一個(gè)數(shù)c1是極坐標(biāo)（1， PI/2），第二個(gè)數(shù)c2是笛卡爾坐標(biāo)（3， 4i），你可以不用考慮轉(zhuǎn)換問題直接運(yùn)算c1（+，-，*，/）c2，運(yùn)算結(jié)果的坐標(biāo)類型由第一個(gè)操作數(shù)決定。

好了，看程序吧。 需要指出的是，為了方便，我把復(fù)數(shù)類的定義和實(shí)現(xiàn)寫在一個(gè)文件里了。

?

?

//

// FILE: my_complex.h

//

?

#ifndef _MY_COMPLEX_H_
#define _MY_COMPLEX_H_

#include <cmath>
#include <cassert>
#include <iostream>

static const double PI = 3.1415926536;
static const double epsilon=1e-8;

class Complex
{
?double real;
?double img;
?bool is_polar;

#define IsZero(x) ?( (bool)((x)<epsilon && (x)>(-epsilon)) )

?void check_zero(void);
public:
?Complex() : real(0.0), img(0.0), is_polar(false) {};

?//no explicit here, otherwise you cannot use codes like Complex X=32; etc.
?Complex(double r) : real(r), img(0.0), is_polar(false) { check_zero(); };

?Complex(double r, double i, bool polar=false) : real(r), img(i), is_polar(polar)
?{ if(is_polar) assert(real>(-epsilon)); check_zero(); };

?inline const double get_real() const;
?inline const double get_img() const;
?inline const double get_rho() const;
?inline const double get_theta() const;
?inline const bool IsPolar() const { return is_polar; };

?void to_polar(void);
?void to_ri(void);
?void conj(void) { img = -img; };

?const Complex& operator = (double a);
?const Complex& operator = (const Complex& a);
?const bool operator == (const Complex& a) const;

?const Complex operator -() const;

?const Complex operator +(const Complex& a) const;
?const Complex operator -(const Complex& a) const;
?const Complex operator *(double a) const;
?const Complex operator *(const Complex& a) const;
?const Complex operator /(double a) const;
?const Complex operator /(const Complex& a) const;

?friend std::ostream& operator <<(std::ostream& os, const Complex& a);

?friend const Complex conj(const Complex& a);
};

void Complex::check_zero(void)
{
?if(IsZero(real))
??real = 0.0;
?if(IsZero(img))
??img = 0.0;
}

void Complex::to_polar(void)
{
?if(is_polar==false)
?{
??double temp_real = get_rho();
??img = get_theta();
??real = temp_real;
??is_polar = true;
?}
}

void Complex::to_ri(void)
{
?if(is_polar)
?{
??double temp_real = get_real();
??img = get_img();
??real = temp_real;
??is_polar = false;
?}
}

const double Complex::get_real() const
{
?return (is_polar ? real*::cos(img):real);
}

const double Complex::get_img() const
{
?return (is_polar ? real*::sin(img):img);
}

const double Complex::get_rho() const
{
?return (is_polar ? real:std::sqrt(real*real+img*img));
}

const double Complex::get_theta() const
{
?if(is_polar)
??return img;
?else
??return ::atan2(img, real);
}

inline const Complex& Complex::operator = (double a)
{
?is_polar = false;
?real = a;
?img = 0;

?return *this;
}

inline const Complex& Complex::operator = (const Complex& a)
{
?is_polar = a.IsPolar();
?real = (is_polar ? a.get_rho():a.get_real());
?img = (is_polar ? a.get_theta():a.get_img());

?return *this;
}

inline const bool Complex::operator == (const Complex& a) const
{
?if(is_polar == a.is_polar)
??return (IsZero(real-a.real) && IsZero(img-a.img));
?else
??return (IsZero(get_real()-a.get_real()) && IsZero(get_img()-a.get_img()));
}

inline const Complex Complex::operator -() const
{
?return (is_polar ? Complex(real, PI+img, true):Complex(-real, -img));
}

const Complex Complex::operator +(const Complex& a) const
{
?Complex temp(get_real()+a.get_real(), get_img()+a.get_img());
?if(is_polar)
??temp.to_polar();
?return temp;
}

const Complex Complex::operator -(const Complex& a) const
{
?Complex temp(get_real()-a.get_real(), get_img()-a.get_img());
?if(is_polar)
??temp.to_polar();
?return temp;
}

const Complex Complex::operator *(double a) const
{
?assert(!IsZero(a));

?if(is_polar)
??return Complex(get_rho()*a, get_theta(), true);
?else
??return Complex(get_real()*a, get_img()*a);
}

const Complex Complex::operator *(const Complex& a) const
{
?if(is_polar)
??return Complex(get_rho()*a.get_rho(), get_theta()+a.get_theta(), true);
?else
??return Complex(get_real()*a.get_real()-get_img()*a.get_img(), get_real()*a.get_img()+get_img()*a.get_real());
}

const Complex Complex::operator /(double a) const
{
?assert(!IsZero(a));

?if(is_polar)
??return Complex(get_rho()/a, get_theta(), true);
?else
??return Complex(get_real()/a, get_img()/a);
}

const Complex Complex::operator /(const Complex& a) const
{
?Complex temp(get_rho()/a.get_rho(), get_theta()-a.get_theta(), true);
?if(is_polar==false)
??temp.to_ri();
?return temp;
}

std::ostream& operator <<(std::ostream& os, const Complex& a)
{
?if(a.IsPolar()==false)
??os<<'('<<a.get_real()<<"r,"<<a.get_img()<<"i)";
?else
??os<<'<'<<a.get_rho()<<"h,"<<a.get_theta()<<"t)";
?return os;
}

inline const Complex conj(const Complex& a)
{
?return Complex(a.real, -a.img, a.is_polar);
}

#endif

?

?

?

?

?

?

?

?

?

?

//

// FILE: fft_dit_base2.cpp

//

?

#include "my_complex.h"

using namespace std;

static const int N=32;

static Complex fft_coeff[N];
static Complex ifft_coeff[N];

void setup_table(void)
{
?for(int i=0; i<N; ++i)
?{
??fft_coeff[i] = Complex(1.0, -PI*i*2/N, true);
??ifft_coeff[i] = Complex(1.0, PI*i*2/N, true);
?}
}

void decimation(Complex data[])
{
?int r=0;
?int i=1;
#define swp(x, y)?{Complex temp=(x); (x)=(y); (y)=temp;}

?while(i<N/2)
?{
??r += N/2;
??swp(data[i], data[r]);
??++i;

??for(int m=N>>1; !((r^=m)&m); m>>=1);
??if(r>i)
??{
???swp(data[i], data[r]);
???swp(data[N-1-i], data[N-1-r]);
??}
??++i;
?}
#undef swp
}

void fft_dit(Complex data[], bool ifft=false)
{
?decimation(data);

?Complex *coeff;
?if(ifft)
??coeff = ifft_coeff;
?else
??coeff = fft_coeff;

?int groups=N/2;
?int step = 2;
?while(groups>=1)
?{
??for(int i=0; i<groups; ++i)
??{
???int first = i*step;
???int second = first+step/2;
???int coeff_idx = 0;
???for(int j=0; j<step/2; ++j)
???{
????Complex temp = data[second]*coeff[coeff_idx];

????data[second] = data[first] - temp;
????data[first]? = data[first] + temp;

????if(ifft==true)
????{
?????data[first] = data[first]/2;
?????data[second] = data[second]/2;
????}

????++first;
????++second;
????coeff_idx += groups;
???}
??}
??groups /= 2;
??step *= 2;
?}
}

int main()
{
?setup_table();

?Complex raw_data[N];
?for(int i=0; i<N; ++i)
??raw_data[i] = i+1;

?cout<<"The data to be fft-ed are:"<<endl;
?for(int i=0; i<N; ++i)
??cout<<raw_data[i]<<' ';
?cout<<endl;

?fft_dit(raw_data);
?cout<<"The fft result are:"<<endl;
?for(int i=0; i<N; ++i)
??cout<<raw_data[i]<<' ';
?cout<<endl;

?fft_dit(raw_data, true);
?cout<<"The ifft-ed data are:"<<endl;
?for(int i=0; i<N; ++i)
??cout<<raw_data[i]<<' ';
?cout<<endl;

?system("PAUSE");
?return 0;
}

?

?

?

?

參考資料:

?

1. http://www.jjj.de/fxt/, 極具參考價(jià)值的網(wǎng)站

?

2. 任何一本講信號(hào)與系統(tǒng),或數(shù)字信號(hào)處理的書

?

總結(jié)

以上是生活随笔為你收集整理的复数基2 DIT FFT程序的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。