From ed1a5f0cfcdabb9c841fbba6403a8790c4d5d881 Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Wed, 11 May 2005 02:36:33 +0000 Subject: [PATCH] made easy to use longs for fixed point --- Makefile | 1 + _kiss_fft_guts.h | 34 ++++++++++----- kiss_fft.h | 6 ++- test/Makefile | 9 +++- test/twotonetest.c | 105 +++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 141 insertions(+), 14 deletions(-) create mode 100644 test/twotonetest.c diff --git a/Makefile b/Makefile index 9e48dc4..9d6b80d 100644 --- a/Makefile +++ b/Makefile @@ -5,6 +5,7 @@ TARBALL=kiss_fft_v$(KFVER).tar.gz ZIPFILE=kiss_fft_v$(KFVER).zip testall: + export DATATYPE=long && cd test && make test export DATATYPE=short && cd test && make test export DATATYPE=float && cd test && make test export DATATYPE=double && cd test && make test diff --git a/_kiss_fft_guts.h b/_kiss_fft_guts.h index f295254..0a1b056 100644 --- a/_kiss_fft_guts.h +++ b/_kiss_fft_guts.h @@ -17,7 +17,7 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND and defines typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ #include "kiss_fft.h" - +#include #define MAXFACTORS 32 /* e.g. an fft of length 128 has 4 factors @@ -42,16 +42,28 @@ struct kiss_fft_state{ C_ADDTO( res , a) : res += a * */ #ifdef FIXED_POINT - -#if defined(CHECK_OVERFLOW) -# define CHECK_OVERFLOW_OP(a,op,b) \ - if ( (long)(a) op (long)(b) > 32767 || (long)(a) op (long)(b) < -32768 ) { \ - fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(long)(a) op (long)(b) ); } +#if (FIXED_POINT==32) +# define FRACBITS 31 +# define SAMPPROD long long +#define SAMP_MAX LONG_MAX +#define SAMP_MIN LONG_MIN +#else +# define FRACBITS 15 +# define SAMPPROD long +#define SAMP_MAX SHRT_MAX +#define SAMP_MIN SHRT_MIN #endif -# define smul(a,b) ( (long)(a)*(b) ) -# define sround( x ) (short)( ( (x) + (1<<14) ) >>15 ) +#if defined(CHECK_OVERFLOW) +# define CHECK_OVERFLOW_OP(a,op,b) \ + if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \ + fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) ); } +#endif + + +# define smul(a,b) ( (SAMPPROD)(a)*(b) ) +# define sround( x ) (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS ) # define S_MUL(a,b) sround( smul(a,b) ) @@ -60,7 +72,7 @@ struct kiss_fft_state{ (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) # define DIVSCALAR(x,k) \ - (x) = sround( smul( x, 32767/k ) ) + (x) = sround( smul( x, SAMP_MAX/k ) ) # define C_FIXDIV(c,div) \ do { DIVSCALAR( (c).r , div); \ @@ -119,8 +131,8 @@ static void kf_cexp(kiss_fft_cpx * x,double phase) /* returns e ** (j*phase) */ { #ifdef FIXED_POINT - x->r = (kiss_fft_scalar) (32767 * cos (phase)); - x->i = (kiss_fft_scalar) (32767 * sin (phase)); + x->r = (kiss_fft_scalar) (SAMP_MAX * cos (phase)); + x->i = (kiss_fft_scalar) (SAMP_MAX * sin (phase)); #else x->r = cos (phase); x->i = sin (phase); diff --git a/kiss_fft.h b/kiss_fft.h index 7446330..a208bc2 100644 --- a/kiss_fft.h +++ b/kiss_fft.h @@ -24,7 +24,11 @@ extern "C" { */ #ifdef FIXED_POINT -# define kiss_fft_scalar short +# if (FIXED_POINT == 32) +# define kiss_fft_scalar long +# else +# define kiss_fft_scalar short +# endif #else # ifndef kiss_fft_scalar /* default is float */ diff --git a/test/Makefile b/test/Makefile index 2226fdb..9ca0dd2 100644 --- a/test/Makefile +++ b/test/Makefile @@ -31,6 +31,11 @@ else SELFTESTSRC=selftest.c endif +ifeq "$(DATATYPE)" "long" + TYPEFLAGS=-DFIXED_POINT=32 -Dkiss_fft_scalar=long + SELFTESTSRC=selftest_long.c +endif + ifeq "$(DATATYPE)" "float" # fftw needs to be built with --enable-float to build this lib FFTWLIB=fftw3f @@ -46,9 +51,9 @@ tools: cd ../tools && make all # for x86 pentium+ machines , these flags work well -#CFLAGS=-Wall -O3 -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer -I.. -I../tools $(WARNINGS) +CFLAGS=-Wall -O3 -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer -I.. -I../tools $(WARNINGS) # If the above flags do not work, try the following -CFLAGS=-Wall -O3 -I.. -I../tools $(WARNINGS) +#CFLAGS=-Wall -O3 -I.. -I../tools $(WARNINGS) $(SELFTEST): $(SELFTESTSRC) $(SRCFILES) $(CC) -o $@ $(CFLAGS) $(TYPEFLAGS) -lm $+ diff --git a/test/twotonetest.c b/test/twotonetest.c new file mode 100644 index 0000000..62fd55d --- /dev/null +++ b/test/twotonetest.c @@ -0,0 +1,105 @@ +#include +#include +#include +#include +#include "kiss_fft.h" +#include "kiss_fftr.h" +#include + +#define pcpx(c)\ + fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) ) + +void fe_fft_test( int N, int bin1,int bin2) +{ + // static kiss_fft_cfg cfg = NULL; + static kiss_fftr_cfg cfg = NULL; + static int lastN = 0; + // kiss_fft_cpx *kin = NULL; + static kiss_fft_cpx *kout = NULL; + static kiss_fft_scalar *kin_r = NULL; + static int halfN = 0; + int flipN; + + int i; + double f1 = (bin1*2.0/N)*M_PI; + double f2 = (bin2*2.0/N)*M_PI; + double sigpow=0; + double noisepow=0; + + if (lastN != N) { + // Free previous structures, if allocated. + if (cfg) + free(cfg); + if (kin_r) + free(kin_r); + if (kout) + free(kout); + + halfN = N / 2; + + // Allocate the KISS FFT configuration structure. + // Although the memory for these structures will be freed + // when the DLL unloads, to be squeaky clean it should be + // freed explicitly. + // cfg = kiss_fft_alloc( N , 0, NULL, NULL); + cfg = kiss_fftr_alloc(N , 0, NULL, NULL); + + // kin = malloc(N * sizeof(kiss_fft_cpx)); + // kin_r = malloc(N * sizeof(kiss_fft_cpx)); + kin_r = malloc(N * sizeof(kiss_fft_scalar)); + kout = malloc(N * sizeof(kiss_fft_cpx)); + + + lastN = N; + } + +#if FIXED_POINT==32 + long maxrange = LONG_MAX; +#else + long maxrange = SHRT_MAX; +#endif + // Convert the floating point "in" array to fixed point. + for (i = 0; i < N; i++) + { + kin_r[i] = (maxrange>>1)*cos(f1*i) + + (maxrange>>1)*cos(f2*i); + } + + // Make the kiss FFT call. + kiss_fftr(cfg, kin_r, kout); + + sigpow = 0; + //printf("["); + for (i=0;i < (N/2+1);++i) { + double tmpr = (double)kout[i].r / (double)maxrange; + double tmpi = (double)kout[i].i / (double)maxrange; + double mag2 = tmpr*tmpr + tmpi*tmpi; + if (i!=0 && i!= N/2) + mag2 *= 2; // all bins except DC and Nyquist have symmetric counterparts implied + sigpow += mag2; + + // subtract out the expected result, anything left is noise + if ( i!=bin1 && i != bin2 ) + noisepow += mag2; + //printf("%d%+di,",(int)kout[i].r,(int)kout[i].i); + } + //printf("]\n"); + printf("TEST %d,%d,%d noise @ %fdB\n",N,bin1,bin2,10*log10(noisepow/sigpow +1e-20) ); + + return(0); +} + + +int main(int argc,char ** argv) +{ + int nfft = 16; + int i,j; + for (i=0;i