#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "FFT.h"

/*----------------------------------------------------
C FFT
By Geoffrey ANNEHEIM
02/10/04
------------------------------------------------------
Informations:
Ce module a été réalisé dans le cadre d'un TPE sur
la reconnaissance vocale. Il se base sur le travail
de Murphy McCauley (MurphyMc@Concentric.NET sur sa
source en Visual Basic.
Ce module contient l'essentiel des fonctions pour
la transformation rapide de fourier (FFT ou TFR).
Les sources d'origines en Turbo Pascal sont
téléchargeables sur:
http://www.fullspectrum.com/deeth
-------------------------------------------------------
Utilisation:
Les fonctions utilisables sont:
FourierTransform() exécute Fast Fourier Transform sur
2 tableaux -- un tableau de réels, un tableau 
d'imaginaires. Si vous n'avez pas de nombres imaginaires
, utilisez un simple tableau intialisé avec la valeur 0.
Cette fonction est capable de faire l'inverse de FFT.
FrequencyOfIndex() peut vous informer de la fréquence
actuelle pointant par un index.
CalcFrequency() transforms une simple fréquence.
--------------------------------------------------------
Notes:
Tous les tableaux doivent commencer par l'index 0.
Le nombre d'échantillons doit être une puissance de 2
c'est-à-dire 2^x.
FrequencyOfIndex() et CalcFrequency() n(ont pas été
testées à fond.
UTILISEZ SES ROUTINES A VOS RISQUES ET PERILES.
-------------------------------------------------------*/

unsigned char NumberOfBitsNeeded (long PowerOfTwo)
{
  for (unsigned char i=0; i<=16; i++) {
    if (PowerOfTwo==(1UL<<i)) return i;}
  return 0;
}

bool IsPowerOfTwo (long x)
{
  if (x<2) return false;
  if (!(x&(x-1))) return true;
  return false;
}

long ReverseBits (long Index, unsigned char NumBits)
{
  long Rev = 0;
  for (unsigned char i=0; i<NumBits; i++)
  {
    Rev = (Rev*2)|(Index&1);
    Index = Index/2;
  } 
  
  return Rev;
}

//bool FFTAudio (long NumSamples, short *RealIn, float *RealOut, float *ImagOut)
bool FFTAudio (long NumSamples, float *RealIn, float *RealOut, float *ImagOut)
{
  /* In this case, NumSamples isn't included (since it's always the same),
     and the imaginary components are left out since they have no meaning here.
     
     I've use Signles instead of Doubles pretty much everywhere. I think this
     makes it faster, but due to type conversion, it actually might not. I should
     check, but I haven't.
     
     The imaginary components have no meaning in this application. I just left out
     the parts of the calculation that need the imaginary input values (which is a
     big speed improvement right there), but we still need output array because
     it's used in the calculation. It's static so that it doesn't get reallocated. 
  */ 
  
  unsigned char NumBits;
  long i, j, k, n;
  long BlockSize, BlockEnd;
  float DeltaAngle, DeltaAr;
  float Alpha, Beta;
  float TR, TI, AR, AI;
  
  /*if ((!IsPowerOfTwo (NumSamples)) || NumSamples<2) {
    printf ("Error in procedure Fourier: NumSamples is %ld, which is not a positive integer power of two.\n", NumSamples);
    return false;}*/
  
  NumBits = NumberOfBitsNeeded (NumSamples);

  for (i=0; i<NumSamples; i++)
  {
    j = ReverseBits (i, NumBits); //I saved time here by pre-calculating all these values
    RealOut [j] = RealIn [i];
    ImagOut [j] = 0;
  }
  
  BlockEnd = 1; BlockSize = 2;
  
  while (BlockSize<=NumSamples)
  {
    DeltaAngle = AngleNumerator/BlockSize;
    
    Alpha = sin(0.5*DeltaAngle);
    Alpha = 2*Alpha*Alpha;
    Beta = sin(DeltaAngle);
    
    i = 0;
    while (i<NumSamples)
    {
      AR = 1; AI = 0;
      
      j = i;
      for (n=0; n<BlockEnd; n++)
      {
        k = j+BlockEnd;
        TR = AR*RealOut [k]-AI*ImagOut [k];
        TI = AI*RealOut [k]+AR*ImagOut [k];
        
        RealOut [k] = RealOut [j]-TR;
        ImagOut [k] = ImagOut [j]-TI;
        RealOut [j] = RealOut [j]+TR;
        ImagOut [j] = ImagOut [j]+TI;
        DeltaAr = Alpha*AR+Beta*AI;
        AI = AI-(Alpha*AI-Beta*AR);
        AR = AR-DeltaAr;
        j++;
      }
      
      i += BlockSize;
    }
    
    BlockEnd = BlockSize;
    BlockSize = BlockSize*2;
  }
  
  return true;
}


/*void iDCT (long size, float *data_in, float *data_out)
{
  unsigned long i, j;
  
  memset (data_out, 0, sizeof(float)*size);
  for (i=0; i<size; i++)
  {
 	for (j=0; j<size; j++)
	  data_out [i] += data_in [j]*cos(M_PI*(float)j*(((float)i+0.5)/(float)size));
	data_out [i] *= 2.0/(float)size;
  }
}*/

void iDCT (long size, float *data_in, float *data_out)
{
  float *din, *dout, k, s = float(size);
  unsigned long i, j;
  
  dout = (float *)memset (data_out, 0, sizeof(float)*size);
  for (i=0; i<size; i++)
  {
    k = (((float)i+0.5)/s);
    for (j=0,din = data_in ; j<size; j++)
      *dout += *din++ * cos(M_PI*(float)j*k);
    *dout++ *= 2.0/(float)size;
  }
}
