Add possibility to use FFTW for OfdmGenerator

3 files changed, 106 insertions, 12 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index f6102e9..b7095c8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -40,6 +40,8 @@ else
 SIMD_CFLAGS =
 endif
 
+bin_PROGRAMS =  odr-dabmod
+
 FFT_DIR=$(top_builddir)/lib/kiss_fft129
 FFT_INC=-I$(FFT_DIR) -I$(FFT_DIR)/tools
 FFT_SRC=$(FFT_DIR)/kiss_fft.c $(FFT_DIR)/kiss_fft.h $(FFT_DIR)/tools/kiss_fftr.c $(FFT_DIR)/tools/kiss_fftr.h kiss_fftsimd.c kiss_fftsimd.h
@@ -47,19 +49,23 @@ FFT_FLG=-ffast-math
 
 .PHONY: kiss_fft129 reed-solomon-4.0
 
-bin_PROGRAMS =  odr-dabmod
-
 DabModulator.cpp:   $(FFT_DIR)
 
 BUILT_SOURCES:  $(FFT_DIR)
 
+if USE_KISS_FFT
+FFT_LDADD=
+else
+FFT_LDADD=-lfftw3f -lm
+endif
+
 $(FFT_DIR):
 	if [ ! -e $(FFT_DIR) ]; then \
 		tar xzf $(top_srcdir)/lib/kiss_fft129.tar.gz -C $(top_builddir)/lib; \
 	fi
 
 odr_dabmod_CPPFLAGS = -Wall $(FFT_INC) $(FFT_FLG) $(SIMD_CFLAGS) $(GITVERSION_FLAGS)
-odr_dabmod_LDADD    = $(ZMQ_LIBS)
+odr_dabmod_LDADD    = $(ZMQ_LIBS) $(FFT_LDADD)
 odr_dabmod_SOURCES  = DabMod.cpp \
                       PcDebug.h \
                       porting.c porting.h \
@@ -114,3 +120,4 @@ EXTRA_DIST	=kiss_fftsimd.c kiss_fftsimd.h
 
 clean-local:
 	rm -rf $(FFT_DIR)
+
diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp
index c7e8ebf..90527d4 100644
--- a/src/OfdmGenerator.cpp
+++ b/src/OfdmGenerator.cpp
@@ -1,6 +1,11 @@
 /*
    Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 Her Majesty
    the Queen in Right of Canada (Communications Research Center Canada)
+
+   Copyright (C) 2014
+   Matthias P. Braendli, matthias.braendli@mpb.li
+
+    http://opendigitalradio.org
  */
 /*
    This file is part of ODR-DabMod.
@@ -21,9 +26,15 @@
 
 #include "OfdmGenerator.h"
 #include "PcDebug.h"
-#include "kiss_fftsimd.h"
+#if USE_FFTW
+#  include "fftw3.h"
+#  define FFT_TYPE fftwf_complex
+#else
+#  include "kiss_fftsimd.h"
+#endif
 
 #include <stdio.h>
+#include <string.h>
 #include <stdexcept>
 #include <assert.h>
 #include <complex>
@@ -37,7 +48,11 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols,
     ModCodec(ModFormat(nbSymbols * nbCarriers * sizeof(FFT_TYPE)),
             ModFormat(nbSymbols * spacing * sizeof(FFT_TYPE))),
     myFftPlan(NULL),
+#if USE_FFTW
+    myFftIn(NULL), myFftOut(NULL),
+#else
     myFftBuffer(NULL),
+#endif
     myNbSymbols(nbSymbols),
     myNbCarriers(nbCarriers),
     mySpacing(spacing)
@@ -57,7 +72,8 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols,
         myNegDst = spacing - (nbCarriers / 2);
         myNegSrc = (nbCarriers + 1) / 2;
         myNegSize = nbCarriers / 2;
-    } else {
+    }
+    else {
         myPosDst = (nbCarriers & 1 ? 0 : 1);
         myPosSrc = nbCarriers / 2;
         myPosSize = (nbCarriers + 1) / 2;
@@ -77,8 +93,25 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols,
     PDEBUG("  myZeroDst: %u\n", myZeroDst);
     PDEBUG("  myZeroSize: %u\n", myZeroSize);
 
+#if USE_FFTW
+    const int N = mySpacing; // The size of the FFT
+    myFftIn = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N);
+    myFftOut = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N);
+    myFftPlan = fftwf_plan_dft_1d(N,
+            myFftIn, myFftOut,
+            FFTW_BACKWARD, FFTW_MEASURE);
+
+    if (sizeof(complexf) != sizeof(FFT_TYPE)) {
+        printf("sizeof(complexf) %d\n", sizeof(complexf));
+        printf("sizeof(FFT_TYPE) %d\n", sizeof(FFT_TYPE));
+        throw std::runtime_error(
+                "OfdmGenerator::process complexf size is not FFT_TYPE size!");
+    }
+#else
     myFftPlan = kiss_fft_alloc(mySpacing, 1, NULL, NULL);
     myFftBuffer = (FFT_TYPE*)memalign(16, mySpacing * sizeof(FFT_TYPE));
+#endif
+
 }
 
 
@@ -86,16 +119,32 @@ OfdmGenerator::~OfdmGenerator()
 {
     PDEBUG("OfdmGenerator::~OfdmGenerator() @ %p\n", this);
 
+#if USE_FFTW
+    if (myFftIn) {
+         fftwf_free(myFftIn);
+    }
+
+    if (myFftOut) {
+         fftwf_free(myFftOut);
+    }
+
+    if (myFftPlan) {
+        fftwf_destroy_plan(myFftPlan);
+    }
+
+#else
     if (myFftPlan != NULL) {
         kiss_fft_free(myFftPlan);
     }
+
     if (myFftBuffer != NULL) {
         free(myFftBuffer);
     }
+
     kiss_fft_cleanup();
+#endif
 }
 
-
 int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut)
 {
     PDEBUG("OfdmGenerator::process(dataIn: %p, dataOut: %p)\n",
@@ -105,6 +154,7 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut)
 
     FFT_TYPE* in = reinterpret_cast<FFT_TYPE*>(dataIn->getData());
     FFT_TYPE* out = reinterpret_cast<FFT_TYPE*>(dataOut->getData());
+
     size_t sizeIn = dataIn->getLength() / sizeof(complexf);
     size_t sizeOut = dataOut->getLength() / sizeof(complexf);
 
@@ -125,7 +175,27 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut)
                 "OfdmGenerator::process output size not valid!");
     }
 
-#ifdef USE_SIMD
+#if USE_FFTW
+    // No SIMD/no-SIMD distinction, it's too early to optimize anything
+    for (size_t i = 0; i < myNbSymbols; ++i) {
+        myFftIn[0][0] = 0;
+        myFftIn[0][1] = 0;
+
+        bzero(&myFftIn[myZeroDst], myZeroSize * sizeof(FFT_TYPE));
+        memcpy(&myFftIn[myPosDst], &in[myPosSrc],
+                myPosSize * sizeof(FFT_TYPE));
+        memcpy(&myFftIn[myNegDst], &in[myNegSrc],
+                myNegSize * sizeof(FFT_TYPE));
+
+        fftwf_execute(myFftPlan);
+
+        memcpy(out, myFftOut, mySpacing * sizeof(FFT_TYPE));
+
+        in += myNbCarriers;
+        out += mySpacing;
+    }
+#else
+#  ifdef USE_SIMD
     for (size_t i = 0, j = 0; i < sizeIn; ) {
         // Pack 4 fft operations
         typedef struct {
@@ -154,7 +224,8 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut)
                     dataBuffer[myNegDst + l].i[k] = cplxIn[i + myNegSrc + l].imag();
                 }
                 i += myNbCarriers;
-            } else {
+            }
+            else {
                 for (size_t l = 0; l < myNbCarriers; ++l) {
                     dataBuffer[l].r[k] = 0.0f;
                     dataBuffer[l].i[k] = 0.0f;
@@ -174,7 +245,7 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut)
             }
         }
     }
-#else
+#  else
     for (size_t i = 0; i < myNbSymbols; ++i) {
         FFT_REAL(myFftBuffer[0]) = 0;
         FFT_IMAG(myFftBuffer[0]) = 0;
@@ -189,7 +260,9 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut)
         in += myNbCarriers;
         out += mySpacing;
     }
+#  endif
 #endif
 
     return sizeOut;
 }
+
diff --git a/src/OfdmGenerator.h b/src/OfdmGenerator.h
index 5be2d03..ec9f14a 100644
--- a/src/OfdmGenerator.h
+++ b/src/OfdmGenerator.h
@@ -1,6 +1,11 @@
 /*
    Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 Her Majesty
    the Queen in Right of Canada (Communications Research Center Canada)
+
+   Copyright (C) 2014
+   Matthias P. Braendli, matthias.braendli@mpb.li
+
+    http://opendigitalradio.org
  */
 /*
    This file is part of ODR-DabMod.
@@ -28,10 +33,14 @@
 
 #include "porting.h"
 #include "ModCodec.h"
-#include "kiss_fftsimd.h"
 
+#if USE_FFTW
+#  include "fftw3.h"
+#else
+#  include "kiss_fftsimd.h"
+#  include <kiss_fft.h>
+#endif
 
-#include <kiss_fft.h>
 #include <sys/types.h>
 
 
@@ -48,8 +57,13 @@ public:
     const char* name() { return "OfdmGenerator"; }
 
 protected:
+#if USE_FFTW
+    fftwf_plan myFftPlan;
+    fftwf_complex *myFftIn, *myFftOut;
+#else
     FFT_PLAN myFftPlan;
     FFT_TYPE *myFftBuffer;
+#endif
     size_t myNbSymbols;
     size_t myNbCarriers;
     size_t mySpacing;
@@ -63,5 +77,5 @@ protected:
     unsigned myZeroSize;
 };
 
-
 #endif // OFDM_GENERATOR_H
+