diff --git a/kiss_fft.c b/kiss_fft.c
index 9fb793e..121c872 100644
--- a/kiss_fft.c
+++ b/kiss_fft.c
@@ -18,6 +18,12 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  fixed or floating point complex numbers.  It also delares the kf_ internal functions.
  */
 
+/* the max length of the sequence of factors of nfft , 
+ MAXFACTORS==32 this should be good up to more 1e15 samples
+ ( 2 * 3**31 )
+ */
+#define MAXFACTORS 32
+
 static void kf_bfly2(
         kiss_fft_cpx * Fout,
         int fstride,
@@ -303,7 +309,7 @@ void * kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
     size_t memneeded = sizeof(kiss_fft_state)
         + sizeof(kiss_fft_cpx)*nfft /* twiddle factors*/
         + sizeof(kiss_fft_cpx)*nfft /* tmpbuf*/
-        + sizeof(int)*nfft /* factors*/
+        + sizeof(int)*2*MAXFACTORS /* factors*/
         + sizeof(kiss_fft_cpx)*nfft; /* scratch*/
 
     if ( lenmem==NULL ) {
diff --git a/test/Makefile b/test/Makefile
index 1c5c030..c0d8ed1 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -13,6 +13,7 @@ BENCHKISS=bm_kiss_$(DATATYPE)
 BENCHFFTW=bm_fftw_$(DATATYPE)
 SELFTEST=st_$(DATATYPE)
 TESTREAL=tr_$(DATATYPE)
+FFTUTIL=kf_$(DATATYPE)
 
 ifeq  "$(DATATYPE)" "short"
  TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
@@ -30,12 +31,15 @@ else
 endif
 
 
-all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL)
+all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL) $(FFTUTIL)
 
 #CFLAGS=-Wall -O3 -ansi -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer
 # If the above flags do not work, try the following
 CFLAGS=-Wall -O3
 
+$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fft2d.c kiss_fftr.c
+	$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
+
 $(SELFTEST): ../kiss_fft.c $(SELFTESTSRC)
 	$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
 
@@ -68,4 +72,4 @@ selftest_short.c:
 	./mk_test.py -s 10 12 14 > selftest_short.c
 
 clean:
-	rm -f *~ bm_* st_* tr_*
+	rm -f *~ bm_* st_* tr_* kf_*
diff --git a/test/compfft.py b/test/compfft.py
new file mode 100755
index 0000000..2811795
--- /dev/null
+++ b/test/compfft.py
@@ -0,0 +1,83 @@
+#!/usr/local/bin/python2.3
+
+# use FFTPACK as a baseline
+import FFT
+from Numeric import *
+import math
+import random
+import sys
+import struct
+
+pi=math.pi
+e=math.e
+j=complex(0,1)
+lims=(-32768,32767)
+
+def randbuf(n,cpx=1):
+    res = array( [ random.uniform( lims[0],lims[1] ) for i in range(n) ] )
+    if cpx:
+        res = res + j*randbuf(n,0)
+    return res
+
+def main():
+    from getopt import getopt
+    import popen2
+    opts,args = getopt( sys.argv[1:],'u:n:Rt:' )
+    opts=dict(opts)
+
+    util = opts.get('-u','./kf_float')
+
+    try:
+        n = int(opts['-n'])
+        cpx = opts.get('-R') is None
+        fmt=opts.get('-t','f')
+    except KeyError:
+        sys.stderr.write("""
+        usage: compfft.py 
+        -n nfft
+        -u utilname : see sample_code/fftutil.c
+        -R : real-optimized version
+        """)
+
+    x = randbuf(n,cpx)
+
+    cmd = '%s -n %d ' % ( util, n )
+    if cpx:
+        xout = FFT.fft(x)
+    else:        
+        cmd += '-R '
+        xout = FFT.real_fft(x)
+
+    proc = popen2.Popen3( cmd , bufsize=len(x) )
+
+    proc.tochild.write( dopack( x , fmt ,cpx ) )
+    proc.tochild.close()
+    xoutcomp = dounpack( proc.fromchild.read( ) , fmt ,1 )
+
+    sigpow = sum( xout * conjugate(xout) )
+    print len(xout)
+    print len(xoutcomp)
+    
+    diff = xout-xoutcomp
+    noisepow = sum( diff * conjugate(diff) )
+
+    snr = 10 * math.log10(abs( sigpow / noisepow ) )
+    print 'NFFT=%d,SNR = %f dB' % (n,snr)
+
+def dopack(x,fmt,cpx):
+    if cpx:
+        s = ''.join( [ struct.pack('ff',c.real,c.imag) for c in x ] )
+    else:
+        s = ''.join( [ struct.pack('f',c) for c in x ] )
+    return s 
+
+def dounpack(x,fmt,cpx):
+    uf = fmt * ( len(x) / 4 )
+    s = struct.unpack(uf,x)
+    if cpx:
+        return array(s[::2]) + array( s[1::2] )*j
+    else:    
+        return array(s )
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/Makefile b/tools/Makefile
index 1c5c030..c0d8ed1 100644
--- a/tools/Makefile
+++ b/tools/Makefile
@@ -13,6 +13,7 @@ BENCHKISS=bm_kiss_$(DATATYPE)
 BENCHFFTW=bm_fftw_$(DATATYPE)
 SELFTEST=st_$(DATATYPE)
 TESTREAL=tr_$(DATATYPE)
+FFTUTIL=kf_$(DATATYPE)
 
 ifeq  "$(DATATYPE)" "short"
  TYPEFLAGS=-DFIXED_POINT -Dkiss_fft_scalar=short
@@ -30,12 +31,15 @@ else
 endif
 
 
-all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL)
+all: $(BENCHKISS) $(SELFTEST) $(BENCHFFTW) $(TESTREAL) $(FFTUTIL)
 
 #CFLAGS=-Wall -O3 -ansi -pedantic -march=pentiumpro -ffast-math -fomit-frame-pointer
 # If the above flags do not work, try the following
 CFLAGS=-Wall -O3
 
+$(FFTUTIL): ../kiss_fft.c fftutil.c kiss_fft2d.c kiss_fftr.c
+	$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
+
 $(SELFTEST): ../kiss_fft.c $(SELFTESTSRC)
 	$(CC) -o $@ $(CFLAGS) -I.. $(TYPEFLAGS) -lm $+
 
@@ -68,4 +72,4 @@ selftest_short.c:
 	./mk_test.py -s 10 12 14 > selftest_short.c
 
 clean:
-	rm -f *~ bm_* st_* tr_*
+	rm -f *~ bm_* st_* tr_* kf_*
diff --git a/tools/fftutil.c b/tools/fftutil.c
index 42d0e69..beb1af5 100644
--- a/tools/fftutil.c
+++ b/tools/fftutil.c
@@ -19,61 +19,90 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 #include <unistd.h>
 
 #include "kiss_fft.h"
+#include "kiss_fft2d.h"
+#include "kiss_fftr.h"
 
-void fft_file(FILE * fin,FILE * fout,int nfft,int nrows,int isinverse,int useascii,int times)
+void fft_file(FILE * fin,FILE * fout,int nfft,int isinverse)
 {
-    int i;
     void *st;
     kiss_fft_cpx * buf;
     kiss_fft_cpx * bufout;
-            
 
-    buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft *nrows );
-    bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft *nrows);
-    if (nrows!=1)
-        st = kiss_fft2d_alloc( nrows,nfft ,isinverse ,0,0);
-    else
-        st = kiss_fft_alloc( nfft ,isinverse ,0,0);
+    buf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
+    bufout = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * nfft );
+    st = kiss_fft_alloc( nfft ,isinverse ,0,0);
 
-    while ( fread( buf , sizeof(kiss_fft_cpx) * nfft * nrows ,1, fin ) > 0 ) {
-        for (i=0;i<times;++i)
-            if (nrows!=1)
-                kiss_fft2d( st , buf ,bufout);
-            else
-                kiss_fft( st , buf ,bufout);
-
-        if (useascii) {
-            int i;
-            for (i=0;i<nfft*nrows;++i) 
-                fprintf(fout, "(%g,%g) ", (double)bufout[i].r,(double)bufout[i].i);
-        }else{
-            fwrite( bufout , sizeof(kiss_fft_cpx) , nfft*nrows , fout );
-        }
+    while ( fread( buf , sizeof(kiss_fft_cpx) * nfft ,1, fin ) > 0 ) {
+        kiss_fft( st , buf ,bufout);
+        fwrite( bufout , sizeof(kiss_fft_cpx) , nfft , fout );
     }
     free(st);
     free(buf);
     free(bufout);
 }
 
+void fft_file2d(FILE * fin,FILE * fout,int nfft,int nrows,int isinverse)
+{
+    void *st;
+    kiss_fft_cpx *buf;
+    kiss_fft_cpx *bufout;
+
+    buf = (kiss_fft_cpx *) malloc (sizeof (kiss_fft_cpx) * nfft * nrows);
+    bufout = (kiss_fft_cpx *) malloc (sizeof (kiss_fft_cpx) * nfft * nrows);
+    st = kiss_fft2d_alloc (nrows, nfft, isinverse, 0, 0);
+
+    while (fread (buf, sizeof (kiss_fft_cpx) * nfft * nrows, 1, fin) > 0) {
+        kiss_fft2d (st, buf, bufout);
+        fwrite (bufout, sizeof (kiss_fft_cpx), nfft * nrows, fout);
+    }
+    free (st);
+    free (buf);
+    free (bufout);
+}
+
+void fft_file_real(FILE * fin,FILE * fout,int nfft,int isinverse)
+{
+    void *st;
+    kiss_fft_scalar * rbuf;
+    kiss_fft_cpx * cbuf;
+
+    rbuf = (kiss_fft_scalar*)malloc(sizeof(kiss_fft_scalar) * nfft );
+    cbuf = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (nfft/2+1) );
+    st = kiss_fftr_alloc( nfft ,isinverse ,0,0);
+
+    if (isinverse==0) {
+        while ( fread( rbuf , sizeof(kiss_fft_scalar) * nfft ,1, fin ) > 0 ) {
+            kiss_fftr( st , rbuf ,cbuf);
+            fwrite( cbuf , sizeof(kiss_fft_cpx) , (nfft/2 + 1) , fout );
+        }
+    }else{
+        while ( fread( cbuf , sizeof(kiss_fft_cpx) * (nfft/2+1) ,1, fin ) > 0 ) {
+            kiss_fftri( st , cbuf ,rbuf);
+            fwrite( rbuf , sizeof(kiss_fft_scalar) , nfft , fout );
+        }
+    }
+    free(st);
+    free(rbuf);
+    free(cbuf);
+}
+
 int main(int argc,char ** argv)
 {
     int nfft=1024;
     int isinverse=0;
+    int isreal=0;
     FILE *fin=stdin;
     FILE *fout=stdout;
-    int useascii=0;
-    int times=1;
     int nrows=1;
 
     while (1) {
-        int c=getopt(argc,argv,"n:iax:r:");
+        int c=getopt(argc,argv,"n:ir:R");
         if (c==-1) break;
         switch (c) {
-            case 'a':useascii=1;break;
             case 'n':nfft = atoi(optarg);break;
             case 'r':nrows = atoi(optarg);break;
             case 'i':isinverse=1;break;
-            case 'x':times=atoi(optarg);break;
+            case 'R':isreal=1;break;
         }
     }
 
@@ -89,7 +118,12 @@ int main(int argc,char ** argv)
         ++optind;
     }
 
-    fft_file(fin,fout,nfft,nrows,isinverse,useascii,times);
+    if (nrows>1 && !isreal)
+        fft_file2d(fin,fout,nfft,nrows,isinverse);
+    else if (nrows==1 && !isreal)
+        fft_file(fin,fout,nfft,isinverse);
+    else if (nrows==1 && isreal)
+        fft_file_real(fin,fout,nfft,isinverse);
 
     if (fout!=stdout) fclose(fout);
     if (fin!=stdin) fclose(fin);