/*
Forward and inverse Discrete Fourier Transform.
Compilation: cc -Wall -Wextra -std=c99 -o dft dft.c -lm

Written by Paweł Sabat.
Web site: sabat.e90.biz

Licensed under Creative Commons Zero license.
See https://creativecommons.org/publicdomain/zero/1.0/ for full license text.
*/

#include <stdio.h>
#include <math.h>

// C99 standard does not define constant PI
#define PI 3.1415926218032837

// structure to represent complex numbers
typedef struct
{
    float r; // real part
    float i; // imaginary part
} cnum;

// complex numbers addition a+b
cnum cnum_add(cnum a, cnum b)
{
    cnum out;

    out.r = a.r + b.r;
    out.i = a.i + b.i;
    return out;
}

// complex numbers multiplication a*b
cnum cnum_mul(cnum a, cnum b)
{
    cnum out;

    out.r = (a.r * b.r) - (a.i * b.i);
    out.i = (a.r * b.i) + (a.i * b.r);
    return out;
}

void dft(cnum *x, /*out*/ cnum *y, int N)
{
    // make sure y is zeroed before use
    for(int n=0; n<N; ++n)
    {
        y[n].r = 0.0;
        y[n].i = 0.0;
    }

    float Nf = (float) N;

    // for complex sinusoid of k-th frequency
    for(int k=0; k<N; ++k)
    {
        // project signal onto that complex sinusoid
        for(int n=0; n<N; ++n)
        {
            cnum csin;
            csin.r = cosf(2*PI*k/Nf*n);
            csin.i = -sinf(2*PI*k/Nf*n); // remember about minus
            cnum addend = cnum_mul(x[n], csin);

            y[k] = cnum_add(y[k], addend);
        }
    }
}

void idft(cnum *y, /*out*/ cnum *x, int N)
{
    // make sure x is zeroed before use
    for(int n=0; n<N; ++n)
    {
        x[n].r = 0.0;
        x[n].i = 0.0;
    }

    float Nf = (float) N;

    // for a given point
    for(int n=0; n<N; ++n)
    {
        // in every complex sinusoid
        for(int k=0; k<N; ++k)
        {
            cnum csin;
            csin.r = cosf(2*PI*k/Nf*n);
            csin.i = sinf(2*PI*k/Nf*n); // no minus this time
            cnum addend = cnum_mul(y[k], csin);

            x[n] = cnum_add(x[n], addend);
        }
    }

    // divide points by N
    for(int n=0; n<N; ++n)
    {
        x[n].r /= Nf;
        x[n].i /= Nf;
    }

}

int main(void)
{
    cnum x[4], y[4], restored_x[4];

    // create some input signal
    x[0].r = 6.0;
    x[0].i = 0.0;
    x[1].r = 3.0;
    x[1].i = 0.0;
    x[2].r = 5.0;
    x[2].i = 0.0;
    x[3].r = 1.0;
    x[3].i = 0.0;

    // compute forward dft
    dft(x, y, 4);

    // display result of forward dft
    printf("%f%+fi\n", y[0].r, y[0].i);
    printf("%f%+fi\n", y[1].r, y[1].i);
    printf("%f%+fi\n", y[2].r, y[2].i);
    printf("%f%+fi\n", y[3].r, y[3].i);

    // compute inverse dft
    idft(y, restored_x, 4);

    // separate result of dft and idft with empty line
    printf("\n");

    // display result of inverse dft
    printf("%f%+fi\n", restored_x[0].r, restored_x[0].i);
    printf("%f%+fi\n", restored_x[1].r, restored_x[1].i);
    printf("%f%+fi\n", restored_x[2].r, restored_x[2].i);
    printf("%f%+fi\n", restored_x[3].r, restored_x[3].i);

    // show the difference between values of original and reconstructed signal
    printf("\n");
    printf("%d %d\n", x[2].r == restored_x[2].r, x[2].i == restored_x[2].i);
    printf("%e, %e\n", x[2].r - restored_x[2].r, x[2].i - restored_x[2].i);

    return 0;
}