/*
CC=================================================================CC
CC                                                                 CC
CC  Subroutine RSRFFT(X,M):                                        CC
CC      A real-valued, in-place, split-radix FFT program           CC
CC      Decimation-in-time, cos/sin in second loop                 CC
CC      and computed recursively                                   CC
CC      Output in order:                                           CC
CC              [ Re(0),Re(1),....,Re(N/2),Im(N/2-1),...Im(1)]     CC
CC                                                                 CC
CC  Input/output                                                   CC
CC      X    Array of input/output (length >= N)                   CC
CC      M    Transform length is N=2**M                            CC
CC                                                                 CC
CC  Calls:                                                         CC
CC      RSTAGE,RBITREV                                             CC
CC                                                                 CC
CC  Author:                                                        CC
CC      H.V. Sorensen,   University of Pennsylvania,  Oct. 1985    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC  Modified:                                                      CC
CC      F. Bonzanigo,    ETH-Zurich,                  Sep. 1986    CC
CC      H.V. Sorensen,   University of Pennsylvania,  Mar. 1987    CC
CC      H.V. Sorensen,   University of Pennsylvania,  Oct. 1987    CC
CC      R. Simard,       Université de Montréal,      Mar. 2014    CC
CC                                                                 CC
CC  Reference:                                                     CC
CC      Sorensen, Jones, Heideman, Burrus :"Real-valued fast       CC
CC      Fourier transform algorithms", IEEE Tran. ASSP,            CC
CC      Vol. ASSP-35, No. 6, pp. 849-864, June 1987                CC
CC                                                                 CC
CC      This program may be used and distributed freely as         CC
CC      as long as this header is included                         CC
CC                                                                 CC
CC=================================================================CC
*/
#include <math.h>

//-------------------Functions Prototypes----------------------------
void rsrfft (double *x, int m);
void rstage (int n, int n2, int n4, double *x1, double *x2, double *x3, double *x4);
void rbitrev (double *x, int m);
//-------------------------------------------------------------------


void rsrfft (double *x, int m)
{
   int n, is, id, i0, n2, n4, k;
   double t1;
   x--;
   n = 1 << m;
   /*C-------Digit reverse counter---------------------------------------C*/
   rbitrev (x, m);
   /*C-----Length two butterflies----------------------------------------C*/
   is = 1;
   id = 4;
   do {
      for (i0 = is; i0 <= n; i0 += id) {
         t1 = x[i0];
         x[i0] = t1 + x[i0 + 1];
         x[i0 + 1] = t1 - x[i0 + 1];
      }
      is = 2 * id - 1;
      id = 4 * id;
   } while (is < n);
   /*C-------L shaped butterflies----------------------------------------C*/
   n2 = 2;
   for (k = 2; k <= m; k++) {
      n2 = n2 * 2;
      n4 = n2 / 4;
      rstage (n, n2, n4, &x[0], &x[n4], &x[2 * n4], &x[3 * n4]);
   }
   return ;
}

/*
C===================================================================C
C  Subroutine RSTAGE - the work-horse of the RFFT                   C
C       Computes a stage of a real-valued split-radix length N      C
C       transform.                                                  C
C  Author                                                           C
C       H.V. Sorensen,   University of Pennsylvania,  Mar. 1987     C
C===================================================================C
*/
void rstage (int n, int n2, int n4, double *x1, double *x2, double *x3, double *x4)
{
   int n8, is, id, i1, i2, j, jn;
   double t1, t2, e, t3, t4, t5;
   double ss1, sd1, ss3, sd3;
   double cc1, cd1, cc3, cd3;
   n8 = n2 / 8;
   is = 0;
   id = n2 * 2;

   do {
      for (i1 = is + 1; i1 <= n; i1 += id) {
         t1 = x4[i1] + x3[i1];
         x4[i1] = x4[i1] - x3[i1];
         x3[i1] = x1[i1] - t1;
         x1[i1] = x1[i1] + t1;
      }
      is = 2 * id - n2;
      id = 4 * id;
   } while (is < n);
   /*C*/
   if (n4 - 1 <= 0)
      return ;

   is = 0;
   id = n2 * 2;
   do {
      for (i2 = is + 1 + n8; i2 <= n; i2 += id) {
         t1 = (x3[i2] + x4[i2]) * .7071067811865475;
         t2 = (x3[i2] - x4[i2]) * .7071067811865475;
         x4[i2] = x2[i2] - t1;
         x3[i2] = -x2[i2] - t1;
         x2[i2] = x1[i2] - t2;
         x1[i2] = x1[i2] + t2;
      }
      is = 2 * id - n2;
      id = 4 * id;
   } while (is < n);
   /*C*/
   if (n8 - 1 <= 0)
      return ;

   e = 2 * 3.14159265358979323 / n2;
   ss1 = sin (e);
   sd1 = ss1;
   sd3 = 3 * sd1 - 4 * (sd1 * sd1 * sd1);
   ss3 = sd3;
   cc1 = cos (e);
   cd1 = cc1;
   cd3 = 4 * (cd1 * cd1 * cd1) - 3. * cd1;
   cc3 = cd3;

   for (j = 2; j <= n8; j++) {
      is = 0;
      id = 2 * n2;
      jn = n4 - 2 * j + 2;
      do {
         for (i1 = is + j; i1 <= n; i1 += id) {
            i2 = i1 + jn;
            t1 = x3[i1] * cc1 + x3[i2] * ss1;
            t2 = x3[i2] * cc1 - x3[i1] * ss1;
            t3 = x4[i1] * cc3 + x4[i2] * ss3;
            t4 = x4[i2] * cc3 - x4[i1] * ss3;
            t5 = t1 + t3;
            t3 = t1 - t3;
            t1 = t2 + t4;
            t4 = t2 - t4;
            x3[i1] = t1 - x2[i2];
            x4[i2] = t1 + x2[i2];
            x3[i2] = -x2[i1] - t3;
            x4[i1] = x2[i1] - t3;
            x2[i2] = x1[i1] - t5;
            x1[i1] = x1[i1] + t5;
            x2[i1] = x1[i2] + t4;
            x1[i2] = x1[i2] - t4;
         }
         is = 2 * id - n2;
         id = 4 * id;
      } while (is < n);
      /*C*/
      t1 = cc1 * cd1 - ss1 * sd1;
      ss1 = cc1 * sd1 + ss1 * cd1;
      cc1 = t1;
      t3 = cc3 * cd3 - ss3 * sd3;
      ss3 = cc3 * sd3 + ss3 * cd3;
      cc3 = t3;
   }
   /*C*/ return  ;
}


/*
CC=================================================================CC
CC                                                                 CC
CC Subroutine RBITREV(X,M):                                        CC
CC      Bitreverses the array X of length 2**M. It generates a     CC
CC      table ITAB (minimum length is SQRT(2**M) if M is even      CC
CC      or SQRT(2*2**M) if M is odd). ITAB need only be generated  CC
CC      once for a given transform length.                         CC
CC                                                                 CC
CC Author:                                                         CC
CC      H.V. Sorensen,   University of Pennsylvania,  Aug. 1987    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC                                                                 CC
CC      This program may be used and distributed freely as long    CC
CC      as this header is included.                                CC
CC                                                                 CC
CC=================================================================CC
*/
void rbitrev (double *x, int m)
{
   int itab[66000];
   int m2, nbit, imax, lbss, i, j, k, l, j0;
   double t1;
   /*C-------Initialization of ITAB array--------------------------------C*/
   m2 = m / 2;
   nbit = 1 << m2;
   if (2 * m2 != m)
      m2 = m2 + 1;
   itab[1] = 0;
   itab[2] = 1;
   imax = 1;

   for (lbss = 2; lbss <= m2; lbss++) {
      imax = 2 * imax;
      for (i = 1; i <= imax; i++) {
         itab[i] = 2 * itab[i];
         itab[i + imax] = 1 + itab[i];
      }
   }

   /*C-------The actual bitreversal--------------------------------------C*/
   for (k = 2; k <= nbit; k++) {
      j0 = nbit * itab[k] + 1;
      i = k;
      j = j0;
      for (l = 2; l <= itab[k] + 1; l++) {
         t1 = x[i];
         x[i] = x[j];
         x[j] = t1;
         i = i + nbit;
         j = j0 + itab[l];
      }
   }
   return ;
}
