// PROGRAM FFTdemo

#include "TrueBASIC.h"
int main();
void coefficients(double g[][2], double fn());
void plotf(double f(), int c);
void FFT(double g[][2], int p);
double fn(double x);
double gr(double x);
double gi(double x);

int N;
double gg[1024][2], eps = 1.0e-13;

int main()
{
    int p, k;
    char Stmp1_[_LBUFF_];
    GWopen(0);

    GWindow(-4.0f*GWaspect(2), -2.0f, 4.0f*GWaspect(2), 2.0f);
    GWsetpen(BLACK, -1, -1, -1);
    GWline((float)(-pi*1.1), 0.0f, (float)(pi*1.1), 0.0f);
    GWline(0.0f, -1.5f, 0.0f, 1.5f);

    plotf(fn, BLACK);

    if(!GWinput("p = ", Stmp1_, _LBUFF_)) {
       GWquit();
       exit(0);
    }
    sscanf(Stmp1_, "%d", &p);

    N = (1 << p);
    coefficients(gg, fn);
    printf("\nN = %d\n", N);
    for(k = -N/2; k <= N/2; ++k) printf("%4d %20.15f %20.15f\n", 
        k, gg[((k+N) % N)][0], gg[((k+N) % N)][1]);
    FFT(gg, p);
    printf("\nN = %d\n", N);
    for(k = -N/2; k <= N/2; ++k) printf("%4d %20.15f %20.15f\n", 
        k, gg[((k+N) % N)][0], gg[((k+N) % N)][1]);
    plotf(gr, RED);
    plotf(gi, BLUE);

    GWquit();
    return 0;
}

void coefficients(double g[][2], double fn())
{
    int k;
    double t, dt;

    t = 0;
    dt = 2*pi/N;
    for(k = 0; k < N; ++k) {
       g[k][0] = fn(t);
       g[k][1] = 0;
       t += dt;
    }
}

void plotf(double f(), int c)
{
    int nplot = 100;
    double t, dt;

    dt = pi/nplot;
    t = -pi;
    GWsetpen(c, -1, -1, -1);
    GWmove2((float)t, (float)(f(t)));
    while(t < pi) {
       t += dt;
       GWline2((float)t, (float)(f(t)));
    }
    GWline2((float)pi, (float)(f(pi)));
}

double fn(double x)
{
#if 0
    return cos(x);
#else
    double f_;
    while(x < 0) x += (2*pi);
    if(fabs(x) < eps)
       f_ = 0;
    else if(x < pi - eps)
       f_ = 1;
    else if(x > pi + eps)
       f_ = -1;
    else
       f_ = 0;
    return f_;
#endif
}

double gr(double x)
{
    double gr_, kx;
    int k;
    while(x < 0) x += (2*pi);
    gr_ = 0;
    for(k = -N/2; k < N/2; ++k) {
       kx = k*x;
       gr_ += (gg[((k+N) % N)][0]*cos(kx)-gg[((k+N) % N)][1]*sin(kx));
    }
    gr_ /= N;
    return gr_;
}

double gi(double x)
{
    double gi_, kx;
    int k;
    while(x < 0) x += (2*pi);
    gi_ = 0;
    for(k = -N/2; k < N/2; ++k) {
       kx = k*x;
       gi_ += (gg[((k+N) % N)][0]*sin(kx)+gg[((k+N) % N)][1]*cos(kx));
    }
    gi_ /= N;
    return gi_;
}