/**************************************************************************
* FFT.c        (creato da Giuseppe Ligorio per Enigma Amiga Run)          *
* Procedure per il calcolo della FFT diretta e inversa                    *
* il programma calcola e stampa un esempio di trasformata di 8 valori     *
* e poi calcola l'antitrasformata degli 8 valori ottenuti per ritornare a *
* quelli originari.                                                       *
**************************************************************************/

/* Inclusione strutture e definizioni di riferimento. */
#include <proto/exec.h>
#include <exec/types.h>
#include <exec/memory.h>
#include <stdio.h>
#include <math.h>

/* Definizione della struttura complex contenenti i dati per un numero complesso */
typedef struct
{
  double real,imag,modulo,angolo;
} complex;

LONG npt = 8;                 /* Numero punti della trasformata */
complex x[] =                 /* Valori da trasformare */
{
  {0.0,0.0,0.0,0.0},          /* Il primo valore non viene considerato perche' */
  {4.0,0.0,0.0,0.0},          /* la trasformazione parte dall'elemento 1 anziche' 0 */
  {2.0,0.0,0.0,0.0},
  {1.0,0.0,0.0,0.0},
  {4.0,0.0,0.0,0.0},
  {6.0,0.0,0.0,0.0},
  {3.0,0.0,0.0,0.0},
  {5.0,0.0,0.0,0.0},
  {2.0,0.0,0.0,0.0}
};

/* Definizione prototipi di funzione. */
void FFT(void);
void InvFFT(void);

/* Trasformazione veloce di Fourier */
void FFT()
{
  register LONG m,irem,l,le,le1,k,ip,i,j;
  register double ur,ui,wr,wi,tr,ti,temp;

  /* Riordinamento dei dati in ingresso */
  j=1;
  for (i=1; i<npt; ++i)
  {
    if (i<j)
    {
      tr = x[j].real; ti = x[j].imag;
      x[j].real = x[i].real;
      x[j].imag = x[i].imag;
      x[i].real = tr; x[i].imag = ti;
    }
    k = npt >> 1;
    while (k<j)
    {
      j -= k;
      k >>= 1;
    }
    j+=k;
  }

 /* Calcolo degli stadi della trasformata */
  m = 0; irem = npt;
  while (irem > 1)
  {
    irem >>= 1;
    m++;
  }
  
  /* Trasformazione */
  for (l=1; l<=m; l++)
  {
    le = (LONG)pow(2.0,(double)l);
    le1 = le >> 1;
    ur = 1.0; ui = 0.0;
    wr = cos(PI/(double)le1);
    wi = -sin(PI/(double)le1);
    for (j=1; j<=le1; ++j)
    {
      i=j;
      while (i<=npt)
      {
        ip = i + le1;
        tr = x[ip].real*ur - x[ip].imag*ui;
        ti = x[ip].imag*ur + x[ip].real*ui;
        x[ip].real = x[i].real - tr;
        x[ip].imag = x[i].imag - ti;
        x[i].real += tr;
        x[i].imag += ti;
        i+=le;
      }
      temp = ur*wr - ui*wi;
      ui = ui*wr + ur*wi;
      ur = temp;
    }
  }
}

/* Trasformazione inversa veloce di Fourier */
void InvFFT()
{
  register LONG m,irem,l,le,le1,k,ip,i,j;
  register double ur,ui,wr,wi,tr,ti,temp;

  j=1;
  for (i=1; i<npt; ++i)
  {
    if (i<j)
    {
      tr = x[j].real; ti = x[j].imag;
      x[j].real = x[i].real;
      x[j].imag = x[i].imag;
      x[i].real = tr; x[i].imag = ti;
    }
    k = npt >> 1;
    while (k<j)
    {
      j -= k;
      k >>= 1;
    }
    j+=k;
  }

  m = 0; irem = npt;
  while (irem > 1)
  {
    irem >>= 1;
    m++;
  }
  
  for (l=1; l<=m; l++)
  {
    le = (LONG)pow(2.0,(double)l);
    le1 = le >> 1;
    ur = 1.0; ui = 0.0;
    wr = cos(PI/(double)le1);
    wi = sin(PI/(double)le1);
    for (j=1; j<=le1; ++j)
    {
      i=j;
      while (i<=npt)
      {
        ip = i + le1;
        tr = x[ip].real*ur - x[ip].imag*ui;
        ti = x[ip].imag*ur + x[ip].real*ui;
        x[ip].real = x[i].real - tr;
        x[ip].imag = x[i].imag - ti;
        x[i].real += tr;
        x[i].imag += ti;
        i+=le;
      }
      temp = ur*wr - ui*wi;
      ui = ui*wr + ur*wi;
      ur = temp;
    }
  }
  for (i=1; i<=npt; ++i)
  {
    x[i].real /= (double)npt;
    x[i].imag /= (double)npt;
  }
}

/* Programma principale. */
void main()
{
  register int kk;
  
  for (kk=1; kk<=npt; kk++)
    printf("%i, %f, %f\n",kk,x[kk].real,x[kk].imag);
  FFT();
  for (kk=1; kk<=npt; kk++)
    printf("%i, %f, %f\n",kk,x[kk].real,x[kk].imag);
  InvFFT();
  for (kk=1; kk<=npt; kk++)
    printf("%i, %f, %f\n",kk,x[kk].real,x[kk].imag);
}