//
//  main.c
//  NaiveFFT
//
//  Created by Jacky Jack on 10/08/2021.
//

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <complex.h>

//data array
#define DATA_SIZE (8)
double data_i[DATA_SIZE] = {1.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};
double data_q[DATA_SIZE] = {0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};
//calc the fft
void fft_if(double *x_i, double *x_q, int n, int inv) {
    
    int i=0,j=0,k=0,m=0, irem=0, sign;
    double tq,ti;
    
    printf("FFT\n");
    k=n/2;
    for (j=1,i=1;i<n; i++) {
        if (i<j) {
            
            //swap values
            ti = x_i[j];
            tq = x_q[j];
            
            x_i[j] = x_i[i];
            x_q[j] = x_q[i];
            
            x_i[i] = ti;
            x_q[i] = tq;
            printf("i=%d j=%d\n",i,j);
            //find middle point
            k = n/2;
            while (k<j) {
                j = j-k;
                k=k/2;
            }
        } else {
            printf("i=%d j=%d\n",i,j);
            k = n/2;
            while (k<j) {
                j = j - k;
                k = k/2;
            }
        }
        j = j+k;
    }
    
    for (i=0;i<DATA_SIZE;i++) {
        printf("%.2f ",x_i[i]);
    }
    printf("\n");
    for (i=0;i<DATA_SIZE;i++) {
        printf("%.2f ",x_q[i]);
    }
    printf("\n");
    
    //calculate number of stages:
    //m=log2(n)
    m=0;
    irem = n;
    while(irem>1) {
        irem = irem/2;
        m = m+1;
    }
    
    
    //FFT or IFFT
    if (inv==1) {
        sign = 1;
    } else {
        sign = -1;
    }
    
    // transform
    for (i=0; i<m; i++) {
        int le=pow(2,i);
        int le1=le/2;
        double ui=1.0;
        double uq=0.0;
        double wi = cos(M_PI/le1);
        double wq = sign*sin(M_PI/le1);
        
        for (j=1; j<le1;j++) {
            i=j;
            while (i<n) {
                int ip = i+le1;
                
                ti = x_i[ip]*ui - x_q[ip]*uq;
                tq = x_q[ip]*ui + x_i[ip]*uq;
                
                x_i[ip] = x_i[i] - ti;
                x_q[ip] = x_q[i] - tq;
                
                x_i[i] = x_i[i] + ti;
                x_q[i] = x_q[i] + tq;
                
                i = i+le;
            }
            double temp = ui*wi - uq*wq;
            uq = uq*wi + ui*wq;
            ui = temp;
        }
    }
    
    //if inverse
    if (inv==1) {
        for (i=0;i<n;i++) {
            x_i[i] = x_i[i]/n;
            x_q[i] = x_q[i]/n;
        }
    }
    
}

//dft works fine
void dft(double *x_i, double *x_q, int n, int inv) {
    double Wn,Wk;
    //static array
    //double Xi[DATA_SIZE],Xq[DATA_SIZE];
    //dynamic array
    double *Xi, *Xq;
    double c,s;
    int i,j;
    
    Xi = malloc(n*sizeof(double));
    Xq = malloc(n*sizeof(double));
    
    Wn = 2*M_PI/n;
    
    if (inv==1) {
        Wn=-Wn;
    }
    
    for (i=0;i<n;i++) {
        Xi[i] = 0.0f;
        Xq[i] = 0.0f;
        
        Wk = i*Wn;
        for (j=0;j<n;j++) {
            c = cos(j*Wk);
            s = sin(j*Wk);
            
            //i - real, q - imaginary
            Xi[i] = Xi[i] + x_i[j]*c + x_q[j]*s;
            Xq[i] = Xq[i] - x_i[j]*s + x_q[j]*c;
        }
        
        if (inv==1) {
            Xi[i] = Xi[i]/n;
            Xq[i] = Xq[i]/n;
        }
    }
    
    for (i=0;i<n;i++) {
        x_i[i] = Xi[i];
        x_q[i] = Xq[i];
    }
    
}

int main(int argc, const char * argv[]) {
    int i;
    
    // insert code here...
    printf("Hello FFT World!\n");
    
    for (i=0;i<DATA_SIZE;i++) {
        printf("%.2f ",data_i[i]);
    }
    printf("\n");
    for (i=0;i<DATA_SIZE;i++) {
        printf("%.2f ",data_q[i]);
    }
    printf("\n");
    fft_if(data_i, data_q, DATA_SIZE, 0);
    
    //dft(data_i, data_q, DATA_SIZE,0);
    
    for (i=0;i<DATA_SIZE;i++) {
        printf("%.2f ",data_i[i]);
    }
    printf("\n");
    for (i=0;i<DATA_SIZE;i++) {
        printf("%.2f ",data_q[i]);
    }
    printf("\n");
    
    return 0;
}