From 206e28f11c5a2745e7045e58025c231f26c51764 Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Wed, 12 Nov 2003 01:09:35 +0000 Subject: [PATCH] made the factorization a separate routine --- fft.py | 39 +++++++++++++++++++++--- kiss_fft.c | 78 ++++++++++++++++++++++++++++++------------------ test/Makefile | 37 ++++++++++------------- test/benchfftw.c | 51 ++++++++++++++++++++++++------- test/pstats.c | 1 - tools/Makefile | 37 ++++++++++------------- 6 files changed, 156 insertions(+), 87 deletions(-) diff --git a/fft.py b/fft.py index a119651..dc82e0b 100644 --- a/fft.py +++ b/fft.py @@ -15,19 +15,19 @@ def fft(f): else: raise Exception('%s not factorable ' % n) - print 'n=%d,p=%d' % (n,p) - print f,' << fin' + #print 'n=%d,p=%d' % (n,p) + #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() diff --git a/kiss_fft.c b/kiss_fft.c index 132ed2a..66b1d3f 100644 --- a/kiss_fft.c +++ b/kiss_fft.c @@ -45,7 +45,7 @@ typedef struct { C_MUL(m,a,b) : m = a*b C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise - C_SUB( res, a,b) : res = a - b + C_SUB( res, a,b) : res = a - b C_SUBFROM( res , a) : res -= a C_ADDTO( res , a) : res += a * */ @@ -86,7 +86,7 @@ typedef struct { kiss_fft_cpx cexp(double phase) /* returns e ** (j*phase) */ { kiss_fft_cpx x; -#ifdef FIXED_POINT +#ifdef FIXED_POINT x.r = (kiss_fft_scalar) ( 32767*cos(phase) ); x.i = (kiss_fft_scalar) ( 32767*sin(phase) ); #else @@ -265,7 +265,7 @@ void bfly5( C_SUB(*Fout1,scratch[5],scratch[6]); C_ADD(*Fout4,scratch[5],scratch[6]); - + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,y2.r) + S_MUL(scratch[8].r,y1.r); scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,y2.r) + S_MUL(scratch[8].i,y1.r); scratch[12].r = - S_MUL(scratch[10].i,y2.i) + S_MUL(scratch[9].i,y1.i); @@ -329,7 +329,7 @@ void fft_work( m=*factors++; for (q=0;q1 && (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) { - /* 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; - } + factor(nfft,st->factors); } /* @@ -426,7 +446,7 @@ void * kiss_fft2d_alloc(int nrows,int ncols,int inverse_fft) st = (kiss_fft2d_state *) malloc ( sizeof(kiss_fft2d_state) + size1 + size2 + sizetmp ); if (!st) return NULL; - + st->minus2 = -2; st->rowst = (kiss_fft_state *)(st+1); /*just beyond kiss_fft2d_state struct */ st->colst = (kiss_fft_state *)( (char*)(st->rowst) + size1 ); @@ -447,7 +467,7 @@ void kiss_fft2d(const void * cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) /*fft each column*/ for (col=0;coltmpbuf[row] = fin[row*ncols + col]; kiss_fft(st->colst,st->tmpbuf); for (row=0;row< nrows ;++row) { @@ -456,7 +476,7 @@ void kiss_fft2d(const void * cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) } /*fft each row */ - for (row=0;row< nrows ;++row) + for (row=0;row< nrows ;++row) kiss_fft(st->rowst , fout + row*ncols ); } diff --git a/test/Makefile b/test/Makefile index c24f3b8..e570fb4 100644 --- a/test/Makefile +++ b/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 - @[ -x ./bm_fftw ] && \ - ./bm_fftw -x $(NUMFFTS) -n $(NFFT) +fftw: $(BENCHFFTW) + ./$(BENCHFFTW) -x $(NUMFFTS) -n $(NFFT) -bm_fftw: benchfftw.c pstats.c - @gcc -o $@ $(CFLAGS) benchfftw.c pstats.c -lm -lfftw3 -L /usr/local/lib/ \ - || echo 'Cannot build FFTW test script' +$(BENCHFFTW): benchfftw.c pstats.c + gcc -o $@ $(CFLAGS) -DDATATYPE$(DATATYPE) benchfftw.c pstats.c -lm -lfftw3f -lfftw3 -L /usr/local/lib/ 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 - for n in $(POW3) ;do \ - ./$(BENCH) -x $(NUMFFTS) -n $$n;\ - [ "$(DATATYPE)" == "double" ] && [ -x ./bm_fftw ] && ./bm_fftw -x $(NUMFFTS) -n $$n || true ; \ +mtime: all $(BENCHFFTW) + for n in $(POW2) $(POW3) $(POW5) ;do \ + echo ============================;\ + ./$(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 diff --git a/test/benchfftw.c b/test/benchfftw.c index 02f79d6..d60b8c7 100644 --- a/test/benchfftw.c +++ b/test/benchfftw.c @@ -4,15 +4,44 @@ #include #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; - fftw_complex * out=NULL; - fftw_plan p; + CPXTYPE * in=NULL; + CPXTYPE * out=NULL; + PLAN p; pstats_init(); @@ -33,28 +62,28 @@ int main(int argc,char ** argv) } } - in=fftw_malloc(sizeof(fftw_complex) * nfft); - out=fftw_malloc(sizeof(fftw_complex) * nfft); + in=FFTMALLOC(sizeof(CPXTYPE) * nfft); + out=FFTMALLOC(sizeof(CPXTYPE) * nfft); for (i=0;i