From 9504aa79c12ab9c4df9e10be4924b1c7c5148a8e Mon Sep 17 00:00:00 2001
From: Mark Borgerding <mark@borgerding.net>
Date: Wed, 15 Oct 2003 02:52:34 +0000
Subject: [PATCH] Fixed generic mixed radix butterfly

'make test' output:
### testing SNR for  1024 point FFTs
#### DOUBLE
snr_t2f = 296.95
snr_f2t = 317.25
#### FLOAT
snr_t2f = 147.96
snr_f2t = 145.14
#### SHORT
snr_t2f = 52.414
snr_f2t = 22.438

#### timing 10000 x 1024 point FFTs
#### DOUBLE
Elapsed:0:03.56 user:2.63 sys:0.19
#### FLOAT
Elapsed:0:01.35 user:1.07 sys:0.10
#### SHORT
Elapsed:0:01.70 user:1.37 sys:0.06
---
 kiss_fft.c | 49 ++++++-------------------------------------------
 1 file changed, 6 insertions(+), 43 deletions(-)

diff --git a/kiss_fft.c b/kiss_fft.c
index 1d5da1d..a8b793d 100644
--- a/kiss_fft.c
+++ b/kiss_fft.c
@@ -37,14 +37,6 @@ typedef struct {
 }kiss_fft_state;
 
 #ifdef FIXED_POINT
-/*
-#define  C_SUB( res, a,b)\
-    do {    (res).r=((a).r-b.r+1)>>1;  (res).i=((a).i-b.i+1)>>1;  }while(0)
-#define C_ADDTO( res , a)\
-    do {    (res).r=((res).r+(a).r+1)>>1;  (res).i=((res).i+(a).i+1)>>1;  }while(0)
-#define C_SUBFROM( res , a)\
-    do {    (res).r=((res).r-(a).r+1)>>1;  (res).i=((res).i-(a).i+1)>>1;  }while(0)
-*/    
     /*  We don't have to worry about overflow from multiplying by twiddle factors since they
      *  all have unity magnitude.  Still need to shift away fractional bits after adding 1/2 for
      *  rounding. */
@@ -84,32 +76,6 @@ void  printcpx(const char * desc,kiss_fft_cpx c)
 {
     printf("%s = (%e,%e)\n",desc,(double)c.r,(double)c.i);
 }
-        
-static
-inline
-kiss_fft_cpx crot(kiss_fft_cpx c, int numquads)
-{
-    kiss_fft_cpx out;
-    switch (numquads&3) {
-        case 0:
-            out.r = c.r;
-            out.i = c.i;
-            break;
-        case 1:
-            out.r = c.i;
-            out.i = -c.r;
-            break;
-        case 2:
-            out.r = -c.r;
-            out.i = -c.i;
-            break;
-        case 3:
-            out.r = -c.i;
-            out.i = c.r;
-            break;
-    }
-    return out;
-}
 
 #define C_ROTADDTO(sum,c,q) \
     do{\
@@ -133,7 +99,6 @@ kiss_fft_cpx crot(kiss_fft_cpx c, int numquads)
         }\
     }while(0) 
 
-
 static 
 void fft_work(
         kiss_fft_cpx * Fout,
@@ -278,13 +243,11 @@ void fft_work(
     kiss_fft_cpx t;
     kiss_fft_cpx * scratch = st->scratch;
     kiss_fft_cpx * twiddles = st->twiddles;
-
-#if 1    
+#if 1
     switch (*factors) {
         case 4: fft4work(Fout,f,fstride,factors,st); return;
         case 2: fft2work(Fout,f,fstride,factors,st); return;
-        default:
-            break;
+        default: break;
     }
 #endif
     p=*factors++;
@@ -295,7 +258,7 @@ void fft_work(
             Fout[q] = *f;
         else
             fft_work( Fout + m*q, f, fstride*p,factors,st);
-        f+= fstride;
+        f += fstride;
     }
 
     for ( u=0; u<m; ++u ) {
@@ -303,15 +266,15 @@ void fft_work(
         for ( q1=0 ; q1<p ; ++q1 ) {
             scratch[q1] = Fout[ k  ];
 #ifdef FIXED_POINT
-            scratch[q1].r >>= 1;
-            scratch[q1].i >>= 1;
+            scratch[q1].r /= p;
+            scratch[q1].i /= p;
 #endif            
             k += m;
         }
 
         k=u;
         for ( q1=0 ; q1<p ; ++q1 ) {
-            int twidx=3;
+            int twidx=0;
             Fout[ k ] = scratch[0];
             for (q=1;q<p;++q ) {
                 int Norig = st->nfft;