/* ----------------------------------------------------------------------   
* Copyright (C) 2010 ARM Limited. All rights reserved.   
*   
* $Date:        15. February 2012  
* $Revision: 	V1.1.0  
*   
* Project: 	    CMSIS DSP Library   
* Title:	    arm_cfft_radix2_f32.c   
*   
* Description:	Radix-2 Decimation in Frequency CFFT & CIFFT Floating point processing function   
*   
*   
* Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
*  
* Version 1.1.0 2012/02/15 
*    Updated with more optimizations, bug fixes and minor API changes.  
*   
* Version 1.0.3 2010/11/29  
*    Initial version   
* -------------------------------------------------------------------- */

#include "arm_math.h"

/**   
 * @ingroup groupTransforms   
 */

/**   
 * @defgroup Radix2_CFFT_CIFFT Radix-2 Complex FFT Functions   
 *   
 * \par   
 * Complex Fast Fourier Transform(CFFT) and Complex Inverse Fast Fourier Transform(CIFFT) is an efficient algorithm to compute Discrete Fourier Transform(DFT) and Inverse Discrete Fourier Transform(IDFT).   
 * Computational complexity of CFFT reduces drastically when compared to DFT.   
 * \par   
 * This set of functions implements CFFT/CIFFT   
 * for Q15, Q31, and floating-point data types.  The functions operates on in-place buffer which uses same buffer for input and output.   
 * Complex input is stored in input buffer in an interleaved fashion.   
 *   
 * \par   
 * The functions operate on blocks of input and output data and each call to the function processes   
 * <code>2*fftLen</code> samples through the transform.  <code>pSrc</code>  points to In-place arrays containing <code>2*fftLen</code> values.   
 * \par  
 * The <code>pSrc</code> points to the array of in-place buffer of size <code>2*fftLen</code> and inputs and outputs are stored in an interleaved fashion as shown below.   
 * <pre> {real[0], imag[0], real[1], imag[1],..} </pre>   
 *   
 * \par Lengths supported by the transform:  
 * \par   
 * Internally, the function utilize a radix-2 decimation in frequency(DIF) algorithm   
 * and the size of the FFT supported are of the lengths [16, 32, 64, 128, 256, 512, 1024, 2048, 4096].  
 *    
 *   
 * \par Algorithm:   
 *   
 * <b>Complex Fast Fourier Transform:</b>   
 * \par    
 * Input real and imaginary data:   
 * <pre>   
 * x(n) = xa + j * ya   
 * x(n+N/2 ) = xb + j * yb   
 * </pre>   
 * where N is length of FFT   
 * \par   
 * Output real and imaginary data:   
 * <pre>   
 * X(2r) = xa'+ j * ya'   
 * X(2r+1) = xb'+ j * yb'   
 * </pre>   
 * \par   
 * Twiddle factors for radix-2 FFT:   
 * <pre>   
 * Wn = cosVal + j * (- sinVal)   
 * </pre>   
 *   
 * \par   
 * \image html CFFT_Radix2.gif "Radix-2 Decimation-in Frequency Complex Fast Fourier Transform"   
 *   
 * \par   
 * Output from Radix-2 CFFT Results in Digit reversal order. Interchange middle two branches of every butterfly results in Bit reversed output.   
 * \par   
 * <b> Butterfly CFFT equations:</b>   
 * <pre>   
 * xa' = xa + xb  
 * ya' = ya + yb  
 * xb' = (xa-xb)* cosVal + (ya-yb) * sinVal   
 * yb' = (ya-yb)* cosVal - (xa-xb) * sinVal   
 * </pre>   
 *   
 *   
 * <b>Complex Inverse Fast Fourier Transform:</b>   
 * \par   
 * CIFFT uses same twiddle factor table as CFFT with modifications in the design equation as shown below.   
 *   
 * \par   
 * <b> Modified Butterfly CIFFT equations:</b>   
 * <pre>   
 * xa' = xa + xb  
 * ya' = ya + yb  
 * xb' = (xa-xb)* cosVal - (ya-yb) * sinVal   
 * yb' = (ya-yb)* cosVal + (xa-xb) * sinVal   
 * </pre>   
 *   
 * \par Instance Structure   
 * A separate instance structure must be defined for each Instance but the twiddle factors and bit reversal tables can be reused.   
 * There are separate instance structure declarations for each of the 3 supported data types.   
 *   
 * \par Initialization Functions   
 * There is also an associated initialization function for each data type.   
 * The initialization function performs the following operations:   
 * - Sets the values of the internal structure fields.   
 * - Initializes twiddle factor table and bit reversal table pointers   
 * \par   
 * Use of the initialization function is optional.   
 * However, if the initialization function is used, then the instance structure cannot be placed into a const data section.   
 * To place an instance structure into a const data section, the instance structure must be manually initialized.   
 * Manually initialize the instance structure as follows:   
 * <pre>   
 *arm_cfft_radix2_instance_f32 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor, onebyfftLen};   
 *arm_cfft_radix2_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};   
 *arm_cfft_radix2_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};   
 * </pre>   
 * \par   
 * where <code>fftLen</code> length of CFFT/CIFFT; <code>ifftFlag</code> Flag for selection of CFFT or CIFFT(Set ifftFlag to calculate CIFFT otherwise calculates CFFT);   
 * <code>bitReverseFlag</code> Flag for selection of output order(Set bitReverseFlag to output in normal order otherwise output in bit reversed order);    
 * <code>pTwiddle</code>points to array of twiddle coefficients; <code>pBitRevTable</code> points to the array of bit reversal table.   
 * <code>twidCoefModifier</code> modifier for twiddle factor table which supports all FFT lengths with same table;    
 * <code>pBitRevTable</code> modifier for bit reversal table which supports all FFT lengths with same table.   
 * <code>onebyfftLen</code> value of 1/fftLen to calculate CIFFT;   
 *  
 * \par Fixed-Point Behavior   
 * Care must be taken when using the fixed-point versions of the CFFT/CIFFT function.   
 * Refer to the function specific documentation below for usage guidelines.   
 */


/**   
 * @addtogroup Radix2_CFFT_CIFFT   
 * @{   
 */

/**   
 * @details   
 * @brief Processing function for the floating-point Radix-2 CFFT/CIFFT.  
 * @param[in]      *S    points to an instance of the floating-point Radix-2 CFFT/CIFFT structure.  
 * @param[in, out] *pSrc points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.  
 * @return none.  
 */

void arm_cfft_radix2_f32(
  const arm_cfft_radix2_instance_f32 * S,
  float32_t * pSrc)
{

  if(S->ifftFlag == 1u)
  {
    /*  Complex IFFT radix-2  */
    arm_radix2_butterfly_inverse_f32(pSrc, S->fftLen, S->pTwiddle,
                                     S->twidCoefModifier, S->onebyfftLen);
  }
  else
  {
    /*  Complex FFT radix-2  */
    arm_radix2_butterfly_f32(pSrc, S->fftLen, S->pTwiddle,
                             S->twidCoefModifier);
  }

  if(S->bitReverseFlag == 1u)
  {
    /*  Bit Reversal */
    arm_bitreversal_f32(pSrc, S->fftLen, S->bitRevFactor, S->pBitRevTable);
  }

}


/**    
 * @} end of Radix2_CFFT_CIFFT group    
 */



/* ----------------------------------------------------------------------    
** Internal helper function used by the FFTs    
** ------------------------------------------------------------------- */

/*    
 * @brief  Core function for the floating-point CFFT butterfly process.   
 * @param[in, out] *pSrc            points to the in-place buffer of floating-point data type.   
 * @param[in]      fftLen           length of the FFT.   
 * @param[in]      *pCoef           points to the twiddle coefficient buffer.   
 * @param[in]      twidCoefModifier twiddle coefficient modifier that supports different size FFTs with the same twiddle factor table.   
 * @return none.   
 */

void arm_radix2_butterfly_f32(
  float32_t * pSrc,
  uint32_t fftLen,
  float32_t * pCoef,
  uint16_t twidCoefModifier)
{

  int i, j, k, l;
  int n1, n2, ia;
  float32_t xt, yt, cosVal, sinVal;

#ifndef ARM_MATH_CM0

  /*  Initializations for the first stage */
  n2 = fftLen;

  n1 = n2;
  n2 = n2 >> 1;
  ia = 0;

  // loop for groups 
  for (i = 0; i < n2; i++)
  {
    cosVal = pCoef[ia * 2];
    sinVal = pCoef[(ia * 2) + 1];

    /*  Twiddle coefficients index modifier */
    ia = ia + twidCoefModifier;

    /*  index calculation for the input as, */
    /*  pSrc[i + 0], pSrc[i + fftLen/1] */
    l = i + n2;

    /*  Butterfly implementation */
    xt = pSrc[2 * i] - pSrc[2 * l];
    pSrc[2 * i] = pSrc[2 * i] + pSrc[2 * l];

    yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
    pSrc[2 * i + 1] = pSrc[2 * l + 1] + pSrc[2 * i + 1];

    pSrc[2u * l] = xt * cosVal + yt * sinVal;

    pSrc[2u * l + 1u] = yt * cosVal - xt * sinVal;

  }                             // groups loop end 

  twidCoefModifier = twidCoefModifier << 1u;

  // loop for stage 
  for (k = fftLen / 2; k > 2; k = k >> 1)
  {
    n1 = n2;
    n2 = n2 >> 1;
    ia = 0;

    // loop for groups 
    for (j = 0; j < n2; j++)
    {
      cosVal = pCoef[ia * 2];
      sinVal = pCoef[(ia * 2) + 1];
      ia = ia + twidCoefModifier;

      // loop for butterfly 
      for (i = j; i < fftLen; i += n1)
      {
        l = i + n2;
        xt = pSrc[2 * i] - pSrc[2 * l];
        pSrc[2 * i] = pSrc[2 * i] + pSrc[2 * l];

        yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
        pSrc[2 * i + 1] = pSrc[2 * l + 1] + pSrc[2 * i + 1];

        pSrc[2u * l] = xt * cosVal + yt * sinVal;

        pSrc[2u * l + 1u] = yt * cosVal - xt * sinVal;

      }                         // butterfly loop end 

    }                           // groups loop end 

    twidCoefModifier = twidCoefModifier << 1u;
  }                             // stages loop end 

  n1 = n2;
  n2 = n2 >> 1;
  ia = 0;

  cosVal = pCoef[ia * 2];
  sinVal = pCoef[(ia * 2) + 1];
  ia = ia + twidCoefModifier;

  // loop for butterfly 
  for (i = 0; i < fftLen; i += n1)
  {
    l = i + n2;
    xt = pSrc[2 * i] - pSrc[2 * l];
    pSrc[2 * i] = (pSrc[2 * i] + pSrc[2 * l]);

    yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
    pSrc[2 * i + 1] = (pSrc[2 * l + 1] + pSrc[2 * i + 1]);

    pSrc[2u * l] = xt;

    pSrc[2u * l + 1u] = yt;

  }                             // groups loop end 

#else

  //N = fftLen; 
  n2 = fftLen;

  // loop for stage 
  for (k = fftLen; k > 1; k = k >> 1)
  {
    n1 = n2;
    n2 = n2 >> 1;
    ia = 0;

    // loop for groups 
    for (j = 0; j < n2; j++)
    {
      cosVal = pCoef[ia * 2];
      sinVal = pCoef[(ia * 2) + 1];
      ia = ia + twidCoefModifier;

      // loop for butterfly 
      for (i = j; i < fftLen; i += n1)
      {
        l = i + n2;
        xt = pSrc[2 * i] - pSrc[2 * l];
        pSrc[2 * i] = pSrc[2 * i] + pSrc[2 * l];

        yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
        pSrc[2 * i + 1] = pSrc[2 * l + 1] + pSrc[2 * i + 1];

        pSrc[2 * l] = (cosVal * xt + sinVal * yt);      // >> 15; 
        pSrc[2 * l + 1] = (cosVal * yt - sinVal * xt);  // >> 15; 

      }
    }
    twidCoefModifier = twidCoefModifier << 1u;
  }

#endif //    #ifndef ARM_MATH_CM0

}


void arm_radix2_butterfly_inverse_f32(
  float32_t * pSrc,
  uint32_t fftLen,
  float32_t * pCoef,
  uint16_t twidCoefModifier,
  float32_t onebyfftLen)
{

  int i, j, k, l;
  int n1, n2, ia;
  float32_t xt, yt, cosVal, sinVal;

#ifndef ARM_MATH_CM0

  //N = fftLen; 
  n2 = fftLen;

  n1 = n2;
  n2 = n2 >> 1;
  ia = 0;

  // loop for groups 
  for (i = 0; i < n2; i++)
  {
    cosVal = pCoef[ia * 2];
    sinVal = pCoef[(ia * 2) + 1];
    ia = ia + twidCoefModifier;

    l = i + n2;
    xt = pSrc[2 * i] - pSrc[2 * l];
    pSrc[2 * i] = pSrc[2 * i] + pSrc[2 * l];

    yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
    pSrc[2 * i + 1] = pSrc[2 * l + 1] + pSrc[2 * i + 1];

    pSrc[2u * l] = xt * cosVal - yt * sinVal;

    pSrc[2u * l + 1u] = yt * cosVal + xt * sinVal;

  }                             // groups loop end 

  twidCoefModifier = twidCoefModifier << 1u;

  // loop for stage 
  for (k = fftLen / 2; k > 2; k = k >> 1)
  {
    n1 = n2;
    n2 = n2 >> 1;
    ia = 0;

    // loop for groups 
    for (j = 0; j < n2; j++)
    {
      cosVal = pCoef[ia * 2];
      sinVal = pCoef[(ia * 2) + 1];
      ia = ia + twidCoefModifier;

      // loop for butterfly 
      for (i = j; i < fftLen; i += n1)
      {
        l = i + n2;
        xt = pSrc[2 * i] - pSrc[2 * l];
        pSrc[2 * i] = pSrc[2 * i] + pSrc[2 * l];

        yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
        pSrc[2 * i + 1] = pSrc[2 * l + 1] + pSrc[2 * i + 1];

        pSrc[2u * l] = xt * cosVal - yt * sinVal;

        pSrc[2u * l + 1u] = yt * cosVal + xt * sinVal;

      }                         // butterfly loop end 

    }                           // groups loop end 

    twidCoefModifier = twidCoefModifier << 1u;
  }                             // stages loop end 

  n1 = n2;
  n2 = n2 >> 1;
  ia = 0;

  cosVal = pCoef[ia * 2];
  sinVal = pCoef[(ia * 2) + 1];
  ia = ia + twidCoefModifier;

  // loop for butterfly 
  for (i = 0; i < fftLen; i += n1)
  {
    l = i + n2;
    xt = pSrc[2 * i] - pSrc[2 * l];
    pSrc[2 * i] = (pSrc[2 * i] + pSrc[2 * l]) * onebyfftLen;

    yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
    pSrc[2 * i + 1] = (pSrc[2 * l + 1] + pSrc[2 * i + 1]) * onebyfftLen;

    pSrc[2u * l] = xt * onebyfftLen;

    pSrc[2u * l + 1u] = yt * onebyfftLen;

  }                             // butterfly loop end 

#else

  //N = fftLen; 
  n2 = fftLen;

  // loop for stage 
  for (k = fftLen; k > 2; k = k >> 1)
  {
    n1 = n2;
    n2 = n2 >> 1;
    ia = 0;

    // loop for groups 
    for (j = 0; j < n2; j++)
    {
      cosVal = pCoef[ia * 2];
      sinVal = pCoef[(ia * 2) + 1];
      ia = ia + twidCoefModifier;

      // loop for butterfly 
      for (i = j; i < fftLen; i += n1)
      {
        l = i + n2;
        xt = pSrc[2 * i] - pSrc[2 * l];
        pSrc[2 * i] = pSrc[2 * i] + pSrc[2 * l];

        yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
        pSrc[2 * i + 1] = pSrc[2 * l + 1] + pSrc[2 * i + 1];

        pSrc[2u * l] = xt * cosVal - yt * sinVal;

        pSrc[2u * l + 1u] = yt * cosVal + xt * sinVal;

      }                         // butterfly loop end 

    }                           // groups loop end 

    twidCoefModifier = twidCoefModifier << 1u;
  }                             // stages loop end 

  n1 = n2;
  n2 = n2 >> 1;
  ia = 0;

  cosVal = pCoef[ia * 2];
  sinVal = pCoef[(ia * 2) + 1];
  ia = ia + twidCoefModifier;

  // loop for butterfly 
  for (i = 0; i < fftLen; i += n1)
  {
    l = i + n2;
    xt = pSrc[2 * i] - pSrc[2 * l];
    pSrc[2 * i] = (pSrc[2 * i] + pSrc[2 * l]) * onebyfftLen;

    yt = pSrc[2 * i + 1] - pSrc[2 * l + 1];
    pSrc[2 * i + 1] = (pSrc[2 * l + 1] + pSrc[2 * i + 1]) * onebyfftLen;

    pSrc[2u * l] = xt * onebyfftLen;

    pSrc[2u * l + 1u] = yt * onebyfftLen;

  }                             // butterfly loop end 

#endif //      #ifndef ARM_MATH_CM0

}