made easy to use longs for fixed point

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

_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 <limits.h>
 
 #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 (FIXED_POINT==32)
-#if defined(CHECK_OVERFLOW)
+# define FRACBITS 31
-#  define CHECK_OVERFLOW_OP(a,op,b)  \
+# define SAMPPROD long long
-	if ( (long)(a) op (long)(b) > 32767 || (long)(a) op (long)(b) < -32768 ) { \
+#define SAMP_MAX LONG_MAX
-		fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(long)(a) op (long)(b) );  }
+#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) )
+#if defined(CHECK_OVERFLOW)
-#   define sround( x )  (short)( ( (x) + (1<<14) ) >>15 )
+#  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->r = (kiss_fft_scalar) (SAMP_MAX * cos (phase));
-    x->i = (kiss_fft_scalar) (32767 * sin (phase));
+    x->i = (kiss_fft_scalar) (SAMP_MAX * sin (phase));
 #else
     x->r = cos (phase);
     x->i = sin (phase);

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 */

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 $+

test/twotonetest.c

@@ -0,0 +1,105 @@
+#include <complex.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include "kiss_fft.h"
+#include "kiss_fftr.h"
+#include <limits.h>
+
+#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<nfft/2;++i) { 
+        for (j=i;j<nfft/2;++j) {
+            fe_fft_test(nfft,i,j);
+        }
+    }
+    printf("sizeof(kiss_fft_scalar) = %d\n",sizeof(kiss_fft_scalar) );
+
+    return 0;
+}