made the factorization a separate routine

fft.py

@@ -15,19 +15,19 @@ def fft(f):
     else:
         raise Exception('%s not factorable ' % n)
 
-    print 'n=%d,p=%d' % (n,p)
+    #print 'n=%d,p=%d' % (n,p)
-    print f,' << fin'
+    #print f,' << fin'
     m = n/p
     Fout=[]
     for q in range(p): # 0,1
         fp = f[q::p]
-        print fp,'<< fp'
+        #print fp,'<< fp'
         Fp = fft( fp )
         Fout.extend( Fp )
 
     for u in range(m):
         scratch = Fout[u::m] # u to end in strides of m
-        print scratch
+        #print scratch
         for q1 in range(p):
             k = q1*m + u  # indices to Fout above that became scratch
             Fout[ k ] = scratch[0] # cuz e**0==1 in loop below
@@ -37,6 +37,37 @@ def fft(f):
 
     return Fout
 
+def real_fft( f ):
+    broken
+
+    N = len(f) / 2
+
+    res = f[::2]
+    ims = f[1::2]
+    fp = [ complex(r,i) for r,i in zip(res,ims) ]
+    Fp = fft( fp )
+
+    Fpr = [ c.real for c in Fp ]
+    Fpi = [ c.imag for c in Fp ]
+
+    F1 = [complex(0,0)]*(N+1)
+    F2 = [complex(0,0)]*(N+1)
+
+    F1[0] = complex( Fpr[0] , 0 )
+    F2[0] = complex( Fpi[0] , 0 )
+    #F1[N] = complex( Fpr[N] , 0 )
+    #F2[N] = complex( Fpi[N] , 0 )
+
+
+    for k in range(1,N):
+        F1[k] = complex( (Fpr[k]+Fpr[N-k])/2 , j*(Fpi[k]-Fpi[N-k])/2 )
+        F2[k] = complex( (Fpi[k]+Fpi[N-k])/2 , j*(Fpr[k]-Fpr[N-k])/2 )
+
+    F = [ complex(0,0) ] * ( N + 1 )
+    for k in range(N+1):
+        F[k] = F1[k] +  e ** ( j*pi*k/N ) * F2[k]
+    return F        
+
 def test(f=range(1024),ntimes=10):
     import time
     t0 = time.time()

kiss_fft.c

@@ -355,9 +355,48 @@ int allocsize(int nfft)
     return allocsize;
 }
 
+/*  factors out powers of 4, powers of 2, then any remaining primes 
+ 
+    facbuf is populated by p1,m1,p2,m2, ...
+
+    where 
+    p[i] * m[i] = m[i-1]
+    m0 = n
+ 
+ *  */
+
+void factor(int n,int * facbuf)
+{
+    int p;
+
+    while ( n>1 && (n&1) == 0) {
+        if ( (n&3) == 0) 
+            p=4;
+        else
+            p=2;
+        n /= p;
+        *facbuf++ = p;
+        *facbuf++ = n;
+    }
+
+    if (n>1) {
+        int floor_sqrt = floor( sqrt( n ) );
+        p=3;
+        do{
+            while (n%p) {
+                p += 2;
+                if ( p>floor_sqrt )
+                    p=n;/* no more factors, skip to end*/
+            }
+            n /= p;
+            *facbuf++ = p;
+            *facbuf++ = n;
+        }while ( n >1);
+    }
+}
+
 void init_state(kiss_fft_state * st,int nfft,int inverse_fft)
 {
-    int nstages=0;
     int i;
     st->nfft=nfft;
     st->inverse = inverse_fft;
@@ -374,26 +413,7 @@ void init_state(kiss_fft_state * st,int nfft,int inverse_fft)
         st->twiddles[i] = cexp( phase );
     }
 
-    while (nfft>1) {
+    factor(nfft,st->factors);
-        /* If you want a new radix, don't forget to put it here */
-        const int divisors[] = {
-		4,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,
-		71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,
-		149,151,157,163,167,173,179,181,191,193,197,199,-1};
-        int p=nfft;
-        i=0;
-        while ( divisors[i] != -1 ) {
-            if ( nfft %  divisors[i] == 0){
-                p =  divisors[i];
-                break;
-            }
-            ++i;
-        }
-        st->factors[2*nstages] = p;
-        nfft /= p;
-        st->factors[2*nstages+1] = nfft;
-        ++nstages;
-    }
 }
 
 /*

test/Makefile

@@ -9,9 +9,10 @@ ifeq "$(DATATYPE)" ""
 endif
 
 UTIL=fftutil_$(DATATYPE)
-BENCH=bm_$(DATATYPE)
+BENCHKISS=bm_kiss_$(DATATYPE)
+BENCHFFTW=bm_fftw_$(DATATYPE)
 
-all: $(UTIL) $(BENCH)
+all: $(UTIL) $(BENCHKISS)
 
 ifeq  "$(DATATYPE)" "short"
  TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
@@ -23,40 +24,34 @@ CFLAGS=-Wall -O3 -ansi -pedantic
 $(UTIL): ../kiss_fft.c fftutil.c 
 	gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) ../kiss_fft.c fftutil.c -lm
 
-$(BENCH): benchkiss.c ../kiss_fft.c pstats.c
+$(BENCHKISS): benchkiss.c ../kiss_fft.c pstats.c
 	gcc -o $@ $(CFLAGS) -I.. benchkiss.c $(TYPEFLAGS) ../kiss_fft.c pstats.c -lm 
 
-fftw: bm_fftw
+fftw: $(BENCHFFTW)
-	@[ -x ./bm_fftw ] && \
+	./$(BENCHFFTW) -x $(NUMFFTS) -n $(NFFT)
-		./bm_fftw -x $(NUMFFTS) -n $(NFFT)
 
-bm_fftw: benchfftw.c pstats.c
+$(BENCHFFTW): benchfftw.c pstats.c
-	@gcc -o $@ $(CFLAGS) benchfftw.c pstats.c -lm -lfftw3 -L /usr/local/lib/ \
+	gcc -o $@ $(CFLAGS) -DDATATYPE$(DATATYPE) benchfftw.c pstats.c -lm -lfftw3f -lfftw3 -L /usr/local/lib/ 
-	 || echo 'Cannot build FFTW test script'
 
 time: all
-	@./$(BENCH) -x $(NUMFFTS) -n $(NFFT) 
+	@./$(BENCHKISS) -x $(NUMFFTS) -n $(NFFT) 
 
-POW2=256 512 1024 2048
+POW2=256 512 1024 2048 4096 8192
 POW3=243 729 2187
 POW5=25 125 625 3125
 
-mtime: all bm_fftw
+mtime: all $(BENCHFFTW)
-	for n in $(POW3) ;do \
+	for n in $(POW2) $(POW3) $(POW5) ;do \
-		./$(BENCH) -x $(NUMFFTS) -n $$n;\
+		echo ============================;\
-		[ "$(DATATYPE)" == "double" ] && [ -x ./bm_fftw ] && ./bm_fftw -x $(NUMFFTS) -n $$n || true ; \
+		./$(BENCHKISS) -x $(NUMFFTS) -n $$n;\
+		[ -x ./$(BENCHFFTW) ] && ./$(BENCHFFTW) -x $(NUMFFTS) -n $$n || true ; \
 	done
 
 snr: all
 	@echo "### testing SNR for $(NFFT) point $(DATATYPE) FFTs"
 	@echo "testkiss( $(NFFT) , '$(DATATYPE)' );" | octave -q
 
-
-ifeq "$(DATATYPE)" "double"
 test: snr time fftw
-else
-test: snr time 
-endif
 
 clean:
-	rm -f *~ fftutil_* bm_*
+	rm -f *~ fftutil_* bm_* *.dat

test/benchfftw.c

@@ -4,15 +4,44 @@
 #include <getopt.h>
 #include "pstats.h"
 
+#ifdef DATATYPEdouble
+
+#define CPXTYPE fftw_complex
+#define PLAN fftw_plan
+#define FFTMALLOC fftw_malloc
+#define MAKEPLAN fftw_plan_dft_1d
+#define DOFFT fftw_execute
+#define DESTROYPLAN fftw_destroy_plan
+#define FFTFREE fftw_free
+
+#elif defined(DATATYPEfloat)
+
+#define CPXTYPE fftwf_complex
+#define PLAN fftwf_plan
+#define FFTMALLOC fftwf_malloc
+#define MAKEPLAN fftwf_plan_dft_1d
+#define DOFFT fftwf_execute
+#define DESTROYPLAN fftwf_destroy_plan
+#define FFTFREE fftwf_free
+
+#endif
+
+#ifndef CPXTYPE
+int main()
+{
+    fprintf(stderr,"Datatype not available in FFTW\n" );
+    return 0;
+}
+#else
 int main(int argc,char ** argv)
 {
     int nfft=1024;
     int isinverse=0;
     int numffts=1000,i;
 
-    fftw_complex * in=NULL;
+    CPXTYPE * in=NULL;
-    fftw_complex * out=NULL;
+    CPXTYPE * out=NULL;
-    fftw_plan p;
+    PLAN p;
 
     pstats_init();
 
@@ -33,28 +62,28 @@ int main(int argc,char ** argv)
       }
     }
 
-    in=fftw_malloc(sizeof(fftw_complex) * nfft);
+    in=FFTMALLOC(sizeof(CPXTYPE) * nfft);
-    out=fftw_malloc(sizeof(fftw_complex) * nfft);
+    out=FFTMALLOC(sizeof(CPXTYPE) * nfft);
     for (i=0;i<nfft;++i ) {
         in[i][0] = rand() - RAND_MAX/2;
         in[i][1] = rand() - RAND_MAX/2;
     }
 
     if ( isinverse )
-        p = fftw_plan_dft_1d(nfft, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
+        p = MAKEPLAN(nfft, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
     else    
-        p = fftw_plan_dft_1d(nfft, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+        p = MAKEPLAN(nfft, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
 
     for (i=0;i<numffts;++i)
-        fftw_execute(p);
+        DOFFT(p);
 
-    fftw_destroy_plan(p);
+    DESTROYPLAN(p);
 
-    fftw_free(in); fftw_free(out);
+    FFTFREE(in); FFTFREE(out);
 
     fprintf(stderr,"fftw\tnfft=%d\tnumffts=%d\n", nfft,numffts);
     pstats_report();
 
     return 0;
 }
-
+#endif

test/pstats.c

@@ -36,6 +36,5 @@ void pstats_report()
 {
     ps_report();
     tms_report();
-    fprintf(stderr,"\n");
 }
 

tools/Makefile

@@ -9,9 +9,10 @@ ifeq "$(DATATYPE)" ""
 endif
 
 UTIL=fftutil_$(DATATYPE)
-BENCH=bm_$(DATATYPE)
+BENCHKISS=bm_kiss_$(DATATYPE)
+BENCHFFTW=bm_fftw_$(DATATYPE)
 
-all: $(UTIL) $(BENCH)
+all: $(UTIL) $(BENCHKISS)
 
 ifeq  "$(DATATYPE)" "short"
  TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
@@ -23,40 +24,34 @@ CFLAGS=-Wall -O3 -ansi -pedantic
 $(UTIL): ../kiss_fft.c fftutil.c 
 	gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) ../kiss_fft.c fftutil.c -lm
 
-$(BENCH): benchkiss.c ../kiss_fft.c pstats.c
+$(BENCHKISS): benchkiss.c ../kiss_fft.c pstats.c
 	gcc -o $@ $(CFLAGS) -I.. benchkiss.c $(TYPEFLAGS) ../kiss_fft.c pstats.c -lm 
 
-fftw: bm_fftw
+fftw: $(BENCHFFTW)
-	@[ -x ./bm_fftw ] && \
+	./$(BENCHFFTW) -x $(NUMFFTS) -n $(NFFT)
-		./bm_fftw -x $(NUMFFTS) -n $(NFFT)
 
-bm_fftw: benchfftw.c pstats.c
+$(BENCHFFTW): benchfftw.c pstats.c
-	@gcc -o $@ $(CFLAGS) benchfftw.c pstats.c -lm -lfftw3 -L /usr/local/lib/ \
+	gcc -o $@ $(CFLAGS) -DDATATYPE$(DATATYPE) benchfftw.c pstats.c -lm -lfftw3f -lfftw3 -L /usr/local/lib/ 
-	 || echo 'Cannot build FFTW test script'
 
 time: all
-	@./$(BENCH) -x $(NUMFFTS) -n $(NFFT) 
+	@./$(BENCHKISS) -x $(NUMFFTS) -n $(NFFT) 
 
-POW2=256 512 1024 2048
+POW2=256 512 1024 2048 4096 8192
 POW3=243 729 2187
 POW5=25 125 625 3125
 
-mtime: all bm_fftw
+mtime: all $(BENCHFFTW)
-	for n in $(POW3) ;do \
+	for n in $(POW2) $(POW3) $(POW5) ;do \
-		./$(BENCH) -x $(NUMFFTS) -n $$n;\
+		echo ============================;\
-		[ "$(DATATYPE)" == "double" ] && [ -x ./bm_fftw ] && ./bm_fftw -x $(NUMFFTS) -n $$n || true ; \
+		./$(BENCHKISS) -x $(NUMFFTS) -n $$n;\
+		[ -x ./$(BENCHFFTW) ] && ./$(BENCHFFTW) -x $(NUMFFTS) -n $$n || true ; \
 	done
 
 snr: all
 	@echo "### testing SNR for $(NFFT) point $(DATATYPE) FFTs"
 	@echo "testkiss( $(NFFT) , '$(DATATYPE)' );" | octave -q
 
-
-ifeq "$(DATATYPE)" "double"
 test: snr time fftw
-else
-test: snr time 
-endif
 
 clean:
-	rm -f *~ fftutil_* bm_*
+	rm -f *~ fftutil_* bm_* *.dat