real fwd and rev fft works

Makefile

@@ -7,5 +7,5 @@ tarball: clean
 
 clean:
 	cd sample_code && make clean
-	rm -f kiss_fft.tar.gz *~ *.pyc kiss_fft.zip fft.py
+	rm -f kiss_fft.tar.gz *~ *.pyc kiss_fft.zip 
 

kiss_fft.h

@@ -59,6 +59,14 @@ void * kiss_fft2d_alloc(int nrows,int ncols,int inverse_fft);
 void kiss_fft2d(const void* cfg_from_alloc , const kiss_fft_cpx *fin,kiss_fft_cpx *fout );
 
 
+/* Real optimized version can save about 40% cpu time vs. complex fft of a real seq.
+ */
+void * kiss_fftr_alloc(int nfft,int inverse_fft);
+void kiss_fftr(const void * cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata);
+void kiss_fftri(const void * cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
+
+
+
 /* when done with the cfg for a given fft size and direction, simply free it*/
 #define kiss_fft_free free
 

kiss_fftr.c

@@ -48,6 +48,8 @@ void * kiss_fftr_alloc(int nfft,int inverse_fft)
 
     for (i=0;i<nfft;++i) {
         double phase = -3.14159265358979323846264338327 * ( (double)i / nfft + .5);
+        if (inverse_fft)
+            phase *= -1;
         st->super_twiddles[i] = kf_cexp( phase );
     }
     return st;
@@ -58,13 +60,13 @@ static void pcpx( kiss_fft_cpx * c)
     printf("%g + %gi\n",c->r,c->i);
 }
 
-void kiss_fftr(const void * cfg,const kiss_fft_scalar *fin,kiss_fft_cpx *fout)
+void kiss_fftr(const void * cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
 {
-    /* input buffer fin is stored row-wise */
+    /* input buffer timedata is stored row-wise */
     kiss_fftr_state *st = ( kiss_fftr_state *)cfg;
     int k,N;
 
-    if (st->minus3 != -3 || st->substate->inverse ) {
+    if ( st->minus3 != -3 || st->substate->inverse) {
         fprintf(stderr,"kiss fft usage error: improper alloc\n");
         exit(1);
     }
@@ -72,10 +74,10 @@ void kiss_fftr(const void * cfg,const kiss_fft_scalar *fin,kiss_fft_cpx *fout)
     N = st->substate->nfft;
 
     /*perform the parallel fft of two real signals packed in real,imag*/
-    kiss_fft( st->substate , (const kiss_fft_cpx*)fin, st->tmpbuf );
+    kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
  
-    fout[0].r = st->tmpbuf[0].r + st->tmpbuf[0].i;
+    freqdata[0].r = st->tmpbuf[0].r + st->tmpbuf[0].i;
-    fout[0].i = 0;
+    freqdata[0].i = 0;
 
     for (k=1;k<N;++k) {
         kiss_fft_cpx fpnk,fpk,f1k;
@@ -86,12 +88,43 @@ void kiss_fftr(const void * cfg,const kiss_fft_scalar *fin,kiss_fft_cpx *fout)
 
         C_ADD( f1k, fpk , fpnk );
         C_SUBFROM( fpk , fpnk );
-        C_MUL( fout[k], fpk , st->super_twiddles[k] );
+        C_MUL( freqdata[k], fpk , st->super_twiddles[k] );
-        C_ADDTO(fout[k],f1k);
+        C_ADDTO(freqdata[k],f1k);
 
-        fout[k].r /= 2;
+        freqdata[k].r /= 2;
-        fout[k].i /= 2;
+        freqdata[k].i /= 2;
     }
+    freqdata[N].r = st->tmpbuf[0].r - st->tmpbuf[0].i;
+    freqdata[N].i = 0;
 }
+    
+void kiss_fftri(const void * cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
+{
+    /* input buffer timedata is stored row-wise */
+    kiss_fftr_state *st = (kiss_fftr_state *) cfg;
+    int k, N;
 
+    if (st->minus3 != -3 || st->substate->inverse == 0) {
+        fprintf (stderr, "kiss fft usage error: improper alloc\n");
+        exit (1);
+    }
 
+    N = st->substate->nfft;
+
+    st->tmpbuf[0].r = freqdata[0].r + freqdata[N].r;
+    st->tmpbuf[0].i = freqdata[0].r - freqdata[N].r;
+    for (k = 1; k <= N / 2; ++k) {
+        kiss_fft_cpx fk, fnkc, fek, fok, tmpbuf;
+        fk = freqdata[k];
+        fnkc.r = freqdata[N - k].r;
+        fnkc.i = -freqdata[N - k].i;
+
+        C_ADD (fek, fk, fnkc);
+        C_SUB (tmpbuf, fk, fnkc);
+        C_MUL (fok, tmpbuf, st->super_twiddles[k]);
+        C_ADD (st->tmpbuf[k], fek, fok);
+        C_SUB (st->tmpbuf[N - k], fek, fok);
+        st->tmpbuf[N - k].i *= -1;
+    }
+    kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
+}

test/test_real.c

@@ -1,24 +1,24 @@
 #include "kiss_fft.h"
-#include <sys/time.h>
+#include <sys/times.h>
+#include <unistd.h>
 
-static double timesnap()
+static double cputime()
 {
-    struct timeval tv;
+    struct tms t;
-    gettimeofday(&tv,NULL);    
+    times(&t);
-    return (double)tv.tv_sec + (double)tv.tv_usec/1000000;
+    return (double)(t.tms_utime + t.tms_stime)/  sysconf(_SC_CLK_TCK) ;
 }
 
-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; 
+    sigpow = noisepow = .00000000000000000001; 
 
     for (k=0;k<n;++k) {
-        sigpow += vec1[k].r * vec1[k].i + 
+        sigpow += vec1[k].r * vec1[k].r + 
                   vec1[k].i * vec1[k].i;
         err = vec1[k].r - vec2[k].r;
         noisepow += err * err;
@@ -36,6 +36,7 @@ double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
 }
 
 #define RANDOM
+ 
 #ifndef RANDOM
 #define NFFT 8
 #else
@@ -43,7 +44,7 @@ double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
 #endif
 
 #ifndef NUMFFTS
-#define NUMFFTS 1000
+#define NUMFFTS 10000
 #endif
 
 void pcpx(const char * msg, kiss_fft_cpx * c)
@@ -60,12 +61,9 @@ int main()
     kiss_fft_cpx cout[NFFT];
     kiss_fft_cpx sout[NFFT];
     
-    const void * kiss_fft_state;
+    void * kiss_fft_state;
-    const void * kiss_fftr_state;
+    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        
@@ -76,27 +74,50 @@ int main()
 /*        printf("in[%d]",i);pcpx("",cin+i); */
     }
 
+    kiss_fft_state = kiss_fft_alloc(NFFT,0);
+    kiss_fftr_state = kiss_fftr_alloc(NFFT,0);
     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) );
+            NFFT,0, snr_compare(cout,sout,(NFFT/2)+1) );
+    free(kiss_fft_state);
+    free(kiss_fftr_state);
 
-    ts = timesnap();
+    kiss_fft_state = kiss_fft_alloc(NFFT,1);
+    kiss_fftr_state = kiss_fftr_alloc(NFFT,1);
+
+    kiss_fft(kiss_fft_state,cout,cin);
+    kiss_fftri(kiss_fftr_state,cout,sin);
+
+    for (i=0;i<NFFT;++i) {
+        sout[i].r = sin[i];
+        sout[i].i = 0;
+        /* printf("sin[%d] = %f\t",i,sin[i]); 
+        printf("cin[%d]",i);pcpx("",cin+i); 
+        printf("sout[%d]",i);pcpx("",sout+i); */ 
+    }
+    
+    printf( "nfft=%d, inverse=%d, snr=%g\n",
+            NFFT,1, snr_compare(cin,cin,NFFT/2) );
+#ifdef RANDOM        
+    ts = cputime();
     for (i=0;i<NUMFFTS;++i) {
         kiss_fft(kiss_fft_state,cin,cout);
     }
-    tfft = timesnap() - ts;
+    tfft = cputime() - ts;
     
-    ts = timesnap();
+    ts = cputime();
     for (i=0;i<NUMFFTS;++i) {
-        kiss_fftr(kiss_fftr_state,sin,cout);
+        /* kiss_fftr(kiss_fftr_state,sin,cout); */
+        kiss_fftri(kiss_fftr_state,cout,sin);
     }
-    trfft = timesnap() - ts;
+    trfft = cputime() - ts;
-    
+
-    printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
+    printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
+#endif        
+    free(kiss_fft_state);
+    free(kiss_fftr_state);
 
-    
     return 0;
 }