#include #include #include #include #include #include #include using namespace std; class Complex{ public: double re, im; Complex(void){ re = im = 0; } inline Complex operator* (double d1){ Complex c; c.re *= d1; c.im *= d1; return c; } inline Complex operator/ (double d0){ Complex c; c.re = re / d0; c.im = im / d0; return c; } inline Complex operator+ (Complex c1){ Complex c; c.re = re + c1.re; c.im = im + c1.im; return c; } inline Complex operator- (Complex c1){ Complex c; c.re = re - c1.re; c.im = im - c1.im; return c; } inline Complex operator* (Complex c1){ Complex c; c.re = re * c1.re - im * c1.im; c.im = re * c1.im + im * c1.re; return c; } inline Complex operator/ (Complex c1){ Complex c; double denomitor = c1.re * c1.re + c1.im * c1.im; c.re = (re * c1.re + im * c1.im) / denomitor; c.im = (- re * c1.im + im * c1.re) / denomitor; return c; } void setPolar(double r, double rad){ re = r * cos(rad); im = r * sin(rad); } double magnitude(void){ return sqrt(re * re + im * im); } }; class Signal{ public: int N, log2N; Complex *W; Signal(int _N){ assert(_N > 2); N = _N; W = new Complex[N]; for(int i = 0; i < N; i++) W[i].setPolar(1, - M_PI * 2 * i / N); } /** DFT */ Complex* dft(Complex *Ft){ Complex *Fk = new Complex[N]; for(int k = 0; k < N; k++){ for(int t = 0; t < N; t++){ Fk[k] = Fk[k] + Ft[t] * W[(k * t) % N]; } } return Fk; } /** IDFT (inverse discrete fourier transform) */ Complex* idft(Complex *Fk){ Complex *Ft = new Complex[N]; for(int t = 0; t < N; t++){ for(int k = 0; k < N; k++) Ft[t] = Ft[t] + Fk[k] * W[N - 1 - (k * t + N - 1) % N]; Ft[t] = Ft[t] / (double)N; } return Ft; } /** FFT wrapper */ Complex* fft(Complex *Ft){ Complex *F = new Complex[N]; for(int i = 0; i < N; i++) F[i] = Ft[i]; fft_core(F, N); int i = 0; for (int j = 1; j < N - 1; j++) { for (int k = N >> 1; k > (i ^= k); k >>= 1); if (i > j) { Complex tmp; tmp = F[j]; F[j] = F[i]; F[i] = tmp; } } return F; } /** FFT core */ void fft_core(Complex *F, int length){ int step = length / 2; int w_step = N / length; for(int i = 0; i < N; i += length){ for(int j = 0; j < step; j++){ Complex c0, c1; c0 = F[i+j]; c1 = F[i+j+step]; F[i+j] = c0 + c1; F[i+j+step] = W[j * w_step] * (c0 - c1); } } if(length > 2) fft_core(F, length/2); } }; void print(Complex *F, int N, string label){ cout << "*** " << label << " ***" << endl; // printf("+++ %s +++\n", label); for(int i = 0; i < N; i++) printf("%5.2f + %5.2f i\n", F[i].re, F[i].im); printf("\n"); } int main(void){ static int N = 8; Complex Ft[N]; for(int i = 0; i < N/4; i++) Ft[i].re = 1, Ft[i].im = 0; for(int i = N/4; i < N/2; i++) Ft[i].re = 0, Ft[i].im = 0; for(int i = N/2; i < N; i++) Ft[i].re = -1, Ft[i].im = 0; print(Ft, N, "Ft"); Signal *sp = new Signal(N); Complex *Fd = sp->dft(Ft); print(Fd, N, "dft"); Complex *Fk = sp->fft(Ft); print(Fk, N, "fft"); Complex *Ft2 = sp->idft(Fk); print(Ft2, N, "idtft"); }