kiss_fftr.c gives correct answers for forward

fft.py

@@ -64,12 +64,15 @@ def real_fft( f,inv ):
     ims = f[1::2]
 
     fp = [ complex(r,i) for r,i in zip(res,ims) ]
+    print 'fft input ', fp
     Fp = fft( fp ,0 )
+    print 'fft output ', Fp
 
     F = []
     for k in range(N):
         F1k = Fp[k] + Fp[-k].conjugate()
         F2k = -j*(Fp[k] - Fp[-k].conjugate())
+        print 'F2k[%d]=%s' % (k,F2k)
 
         F.append( ( F1k + e ** ( -j*pi*k/N ) * F2k ) * .5 )
 

test/Makefile

@@ -11,10 +11,11 @@ ifeq "$(DATATYPE)" ""
 endif
 
 UTIL=fftutil_$(DATATYPE)
+TESTREAL=tr_$(DATATYPE)
 BENCHKISS=bm_kiss_$(DATATYPE)
 BENCHFFTW=bm_fftw_$(DATATYPE)
 
-all: $(UTIL) $(BENCHKISS)
+all: $(UTIL) $(BENCHKISS) $(TESTREAL)
 
 ifeq  "$(DATATYPE)" "short"
  TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
@@ -29,6 +30,9 @@ CFLAGS=-Wall -O3 -ansi -pedantic -march=pentiumpro -ffast-math -fomit-frame-poin
 $(UTIL): ../kiss_fft.c ../kiss_fft2d.c fftutil.c
 	gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) ../kiss_fft.c ../kiss_fft2d.c fftutil.c -lm
 
+$(TESTREAL): ../kiss_fft.c ../kiss_fftr.c test_real.c
+	gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
+
 $(BENCHKISS): benchkiss.c ../kiss_fft.c pstats.c
 	gcc -o $@ $(CFLAGS) -I.. benchkiss.c $(TYPEFLAGS) ../kiss_fft.c pstats.c -lm 
 

test/mk_test.py

@@ -20,13 +20,23 @@ def c_format(v,round=0):
         s= ','.join( [ '{%.60f ,%.60f }' %(c.real,c.imag) for c in v ] ) 
         return re.sub(r'\.?0+ ',' ',s)
 
-def test_vec( n,inverse ):
+def test_cpx( n,inverse ,short):
     v = randvec(n,1)
+    scale = 1
+
     if inverse:
         tvecout = FFT.inverse_fft(v)
-        tvecout = [ c * len(v) for c in tvecout ]
+        if short:
+            scale = 1
+        else:            
+            scale = len(v)
     else:
         tvecout = FFT.fft(v)
+        if short:
+            scale = 1.0/len(v)
+
+    tvecout = [ c * scale for c in tvecout ]
+
 
     s="""#define NFFT %d""" % len(v) + """
     {
@@ -75,7 +85,13 @@ double snr_compare( kiss_fft_cpx * test_vec_out,kiss_fft_cpx * testbuf, int n)
     return s
 
 def main():
-    fftsizes = sys.argv[1:]
+
+    from getopt import getopt
+    opts,args = getopt(sys.argv[1:],'s')
+    opts = dict(opts)
+    short = int( opts.has_key('-s') )
+
+    fftsizes = args
     if not fftsizes:
         fftsizes = [ 1800 ]
     print '#include "kiss_fft.h"'
@@ -83,8 +99,8 @@ def main():
     print "int main() {"
     for n in fftsizes:
         n = int(n)
-        print test_vec(n,0)
+        print test_cpx(n,0,short)
-        print test_vec(n,1)
+        print test_cpx(n,1,short)
     print """
     return 0;
 }

test/test_real.c

@@ -0,0 +1,75 @@
+#include "kiss_fft.h"
+
+
+void * kiss_fftr_alloc(int nfft,int inverse_fft);
+void kiss_fftr(const void * cfg,const kiss_fft_scalar *fin,kiss_fft_cpx *fout);
+
+double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
+{
+    int k;
+    double sigpow,noisepow,err,snr,scale=0;
+    sigpow = noisepow = .000000000000000000000000000001; 
+
+    for (k=0;k<n;++k) {
+        sigpow += vec1[k].r * vec1[k].i + 
+                  vec1[k].i * vec1[k].i;
+        err = vec1[k].r - vec2[k].r;
+        noisepow += err * err;
+        err = vec1[k].i - vec2[k].i;
+        noisepow += err * err;
+
+        if (vec1[k].r)
+            scale += vec2[k].r / vec1[k].r;
+    }
+    snr = 10*log10( sigpow / noisepow );
+    scale /= n;
+    if (snr<10)
+        printf( "\npoor snr, try a scaling factor %f\n" , scale );
+    return snr;
+}
+
+#define RANDOM
+#ifndef RANDOM
+#define NFFT 8
+#else
+#define NFFT 120
+#endif
+
+void pcpx(const char * msg, kiss_fft_cpx * c)
+{
+    printf("%s: %g + %gi\n",msg,c->r,c->i);
+}
+
+int main()
+{
+    int i;
+    kiss_fft_cpx cin[NFFT];
+    kiss_fft_scalar sin[NFFT] = {0.309655,0.815653,0.768570,0.591841,0.404767,0.637617,0.007803,0.012665};
+    kiss_fft_cpx cout[NFFT];
+    kiss_fft_cpx sout[NFFT];
+    
+    const void * kiss_fft_state;
+    const void * kiss_fftr_state;
+    int inverse = 0;
+    
+    kiss_fft_state = kiss_fft_alloc(NFFT,inverse);
+    kiss_fftr_state = kiss_fftr_alloc(NFFT,inverse);
+
+    for (i=0;i<NFFT;++i) {
+#ifdef RANDOM        
+        sin[i] = (kiss_fft_scalar)(rand()-RAND_MAX/2);
+#endif        
+        cin[i].r = sin[i];
+        cin[i].i = 0;
+/*        printf("in[%d]",i);pcpx("",cin+i); */
+    }
+
+    kiss_fft(kiss_fft_state,cin,cout);
+    kiss_fftr(kiss_fftr_state,sin,sout);
+
+    printf( "nfft=%d, inverse=%d, snr=%g\n",
+            NFFT,inverse, snr_compare(cout,sout,NFFT/2) );
+    
+    return 0;
+}
+

tools/Makefile

@@ -11,10 +11,11 @@ ifeq "$(DATATYPE)" ""
 endif
 
 UTIL=fftutil_$(DATATYPE)
+TESTREAL=tr_$(DATATYPE)
 BENCHKISS=bm_kiss_$(DATATYPE)
 BENCHFFTW=bm_fftw_$(DATATYPE)
 
-all: $(UTIL) $(BENCHKISS)
+all: $(UTIL) $(BENCHKISS) $(TESTREAL)
 
 ifeq  "$(DATATYPE)" "short"
  TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
@@ -29,6 +30,9 @@ CFLAGS=-Wall -O3 -ansi -pedantic -march=pentiumpro -ffast-math -fomit-frame-poin
 $(UTIL): ../kiss_fft.c ../kiss_fft2d.c fftutil.c
 	gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) ../kiss_fft.c ../kiss_fft2d.c fftutil.c -lm
 
+$(TESTREAL): ../kiss_fft.c ../kiss_fftr.c test_real.c
+	gcc -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
+
 $(BENCHKISS): benchkiss.c ../kiss_fft.c pstats.c
 	gcc -o $@ $(CFLAGS) -I.. benchkiss.c $(TYPEFLAGS) ../kiss_fft.c pstats.c -lm