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;
}