From 8c5752567cf9acb6d025653a87cc75e713da443f Mon Sep 17 00:00:00 2001 From: Georgi Gerganov Date: Sat, 28 Nov 2020 09:40:46 +0200 Subject: [PATCH] remove FFTW stuff using a simple FFT algorithm instead --- .gitignore | 2 - CMakeLists.txt | 5 +- README.md | 2 +- cmake/FindFFTW.cmake | 27 ------- compile.sh | 2 +- main.cpp | 178 +++++++++++++++++++++++++++++-------------- 6 files changed, 123 insertions(+), 93 deletions(-) delete mode 100644 cmake/FindFFTW.cmake diff --git a/.gitignore b/.gitignore index 358f547..40d145e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,5 @@ build build_timestamp.h compile_* -fftw-3.3.3 *.js *.wasm -lib diff --git a/CMakeLists.txt b/CMakeLists.txt index b3e237e..556f168 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,7 +13,6 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -W -Wall -Wno-long-long -peda # ## Dependencies find_package(Threads REQUIRED) -find_package(FFTW REQUIRED) find_package(SDL2) if (NOT USE_FINDSDL2 AND NOT SDL2_FOUND) @@ -28,5 +27,5 @@ endif() string(STRIP "${SDL2_LIBRARIES}" SDL2_LIBRARIES) add_executable(wave-share main.cpp) -target_include_directories(wave-share PUBLIC ${FFTW_INCLUDE_DIRS} ${SDL2_INCLUDE_DIRS}) -target_link_libraries(wave-share PUBLIC ${CMAKE_THREAD_LIBS_INIT} ${FFTW_LIBRARIES} ${SDL2_LIBRARIES}) +target_include_directories(wave-share PUBLIC ${SDL2_INCLUDE_DIRS}) +target_link_libraries(wave-share PUBLIC ${CMAKE_THREAD_LIBS_INIT} ${SDL2_LIBRARIES}) diff --git a/README.md b/README.md index 1a09f06..22926eb 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,7 @@ For convenience, a [simple WebRTC hack](https://github.com/diafygi/webrtc-ips) i ### Web Assembly module `wave.wasm` -You will need an Emscripten compiler. Additionally, you need [FFTW](http://www.fftw.org) built with Emscripten. Run the ``compile.sh`` script. +You will need an Emscripten compiler. Run the ``compile.sh`` script. ### CLI tool `wave-share` diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake deleted file mode 100644 index 941bec6..0000000 --- a/cmake/FindFFTW.cmake +++ /dev/null @@ -1,27 +0,0 @@ -# - Find FFTW -# Find the native FFTW includes and library -# -# FFTW_INCLUDE_DIRS - where to find fftw3.h -# FFTW_LIBRARIES - List of libraries when using FFTW. -# FFTWF_LIBRARIES - List of libraries when using FFTW single precision. -# FFTW_FOUND - True if FFTW found. - -if (FFTW_INCLUDE_DIRS) - # Already in cache, be silent - set (FFTW_FIND_QUIETLY TRUE) -endif (FFTW_INCLUDE_DIRS) - -find_path (FFTW_INCLUDE_DIRS fftw3.h) - -#find_library (FFTW_LIBRARIES NAMES fftw3f) - -find_library (FFTW_LIBRARIES NAMES fftw3) -find_library (FFTWF_LIBRARIES NAMES fftw3f) -set (FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTWF_LIBRARIES}) - -# handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if -# all listed variables are TRUE -include (FindPackageHandleStandardArgs) -find_package_handle_standard_args (FFTW DEFAULT_MSG FFTW_LIBRARIES FFTW_INCLUDE_DIRS) - -mark_as_advanced (FFTW_LIBRARIES FFTW_INCLUDE_DIRS) diff --git a/compile.sh b/compile.sh index 7c7bfcc..2543dfb 100755 --- a/compile.sh +++ b/compile.sh @@ -2,7 +2,7 @@ echo "static const char * BUILD_TIMESTAMP=\"`date`\";" > build_timestamp.h -em++ -Wall -Wextra -O3 -std=c++11 -s USE_SDL=2 -s WASM=1 ./main.cpp -o wave.js -I ./fftw-3.3.3/api ./lib/libfftw3f.a \ +em++ -Wall -Wextra -O3 -std=c++11 -s USE_SDL=2 -s WASM=1 ./main.cpp -o wave.js \ -s EXPORTED_FUNCTIONS='["_getText", "_getSampleRate", "_setText", "_getAverageRxTime_ms", "_setParameters", "_getFramesLeftToRecord", "_getFramesToRecord", "_getFramesLeftToAnalyze", "_getFramesToAnalyze", diff --git a/main.cpp b/main.cpp index 44ab6bf..2696d6c 100644 --- a/main.cpp +++ b/main.cpp @@ -3,7 +3,6 @@ * \author Georgi Gerganov */ -#include "fftw3.h" #include "reed-solomon/rs.hpp" #include @@ -17,6 +16,7 @@ #include #include #include +#include #ifndef M_PI #define M_PI 3.14159265358979323846f @@ -48,57 +48,126 @@ struct DataRxTx; static DataRxTx *g_data = nullptr; namespace { - //constexpr float IRAND_MAX = 1.0f/RAND_MAX; - //inline float frand() { return ((float)(rand()%RAND_MAX)*IRAND_MAX); } - constexpr double kBaseSampleRate = 48000.0; - constexpr auto kMaxSamplesPerFrame = 1024; - constexpr auto kMaxDataBits = 256; - constexpr auto kMaxDataSize = 256; - constexpr auto kMaxLength = 140; - constexpr auto kMaxSpectrumHistory = 4; - constexpr auto kMaxRecordedFrames = 64*10; - constexpr auto kDefaultFixedLength = 82; +constexpr double kBaseSampleRate = 48000.0; +constexpr auto kMaxSamplesPerFrame = 1024; +constexpr auto kMaxDataBits = 256; +constexpr auto kMaxDataSize = 256; +constexpr auto kMaxLength = 140; +constexpr auto kMaxSpectrumHistory = 4; +constexpr auto kMaxRecordedFrames = 64*10; +constexpr auto kDefaultFixedLength = 82; - enum TxMode { - FixedLength = 0, - VariableLength, - }; +// FFT routines taken from https://stackoverflow.com/a/37729648/4039976 - using AmplitudeData = std::array; - using AmplitudeData16 = std::array; - using SpectrumData = std::array; - using RecordedData = std::array; +int log2(int N) { + int k = N, i = 0; + while(k) { + k >>= 1; + i++; + } + return i - 1; +} - inline void addAmplitudeSmooth(const AmplitudeData & src, AmplitudeData & dst, float scalar, int startId, int finalId, int cycleMod, int nPerCycle) { - int nTotal = nPerCycle*finalId; - float frac = 0.15f; - float ds = frac*nTotal; - float ids = 1.0f/ds; - int nBegin = frac*nTotal; - int nEnd = (1.0f - frac)*nTotal; - for (int i = startId; i < finalId; i++) { - float k = cycleMod*finalId + i; - if (k < nBegin) { - dst[i] += scalar*src[i]*(k*ids); - } else if (k > nEnd) { - dst[i] += scalar*src[i]*(((float)(nTotal) - k)*ids); - } else { - dst[i] += scalar*src[i]; +int check(int n) { + return n > 0 && (n & (n - 1)) == 0; +} + +int reverse(int N, int n) { + int j, p = 0; + for(j = 1; j <= log2(N); j++) { + if(n & (1 << (log2(N) - j))) + p |= 1 << (j - 1); + } + return p; +} + +void ordina(std::complex* f1, int N) { + std::complex f2[kMaxSamplesPerFrame]; + for(int i = 0; i < N; i++) + f2[i] = f1[reverse(N, i)]; + for(int j = 0; j < N; j++) + f1[j] = f2[j]; +} + +void transform(std::complex* f, int N) { + ordina(f, N); //first: reverse order + std::complex *W; + W = (std::complex *)malloc(N / 2 * sizeof(std::complex)); + W[1] = std::polar(1., -2. * M_PI / N); + W[0] = 1; + for(int i = 2; i < N / 2; i++) + W[i] = pow(W[1], i); + int n = 1; + int a = N / 2; + for(int j = 0; j < log2(N); j++) { + for(int i = 0; i < N; i++) { + if(!(i & n)) { + std::complex temp = f[i]; + std::complex Temp = W[(i * a) % (n * a)] * f[i + n]; + f[i] = temp + Temp; + f[i + n] = temp - Temp; } } + n *= 2; + a = a / 2; } + free(W); +} - template - float getTime_ms(const T & tStart, const T & tEnd) { - return ((float)(std::chrono::duration_cast(tEnd - tStart).count()))/1000.0; +void FFT(std::complex* f, int N, float d) { + transform(f, N); + for(int i = 0; i < N; i++) + f[i] *= d; //multiplying by step +} + +void FFT(float * src, std::complex* dst, int N, float d) { + for (int i = 0; i < N; ++i) { + dst[i].real(src[i]); + dst[i].imag(0); + } + FFT(dst, N, d); +} + +enum TxMode { + FixedLength = 0, + VariableLength, +}; + +using AmplitudeData = std::array; +using AmplitudeData16 = std::array; +using SpectrumData = std::array; +using RecordedData = std::array; + +inline void addAmplitudeSmooth(const AmplitudeData & src, AmplitudeData & dst, float scalar, int startId, int finalId, int cycleMod, int nPerCycle) { + int nTotal = nPerCycle*finalId; + float frac = 0.15f; + float ds = frac*nTotal; + float ids = 1.0f/ds; + int nBegin = frac*nTotal; + int nEnd = (1.0f - frac)*nTotal; + for (int i = startId; i < finalId; i++) { + float k = cycleMod*finalId + i; + if (k < nBegin) { + dst[i] += scalar*src[i]*(k*ids); + } else if (k > nEnd) { + dst[i] += scalar*src[i]*(((float)(nTotal) - k)*ids); + } else { + dst[i] += scalar*src[i]; } - - int getECCBytesForLength(int len) { - return std::max(4, 2*(len/5)); } } +template +float getTime_ms(const T & tStart, const T & tEnd) { + return ((float)(std::chrono::duration_cast(tEnd - tStart).count()))/1000.0; +} + +int getECCBytesForLength(int len) { + return std::max(4, 2*(len/5)); +} +} + struct DataRxTx { DataRxTx(int aSampleRateOut, int aSampleRate, int aSamplesPerFrame, int aSampleSizeB, const char * text) { sampleSizeBytes = aSampleSizeB; @@ -213,17 +282,9 @@ struct DataRxTx { rxData.fill(0); - if (fftPlan) fftwf_destroy_plan(fftPlan); - if (fftIn) fftwf_free(fftIn); - if (fftOut) fftwf_free(fftOut); - - fftIn = (float*) fftwf_malloc(sizeof(float)*samplesPerFrame); - fftOut = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*samplesPerFrame); - fftPlan = fftwf_plan_dft_r2c_1d(1*samplesPerFrame, fftIn, fftOut, FFTW_ESTIMATE); - for (int i = 0; i < samplesPerFrame; ++i) { - fftOut[i][0] = 0.0f; - fftOut[i][1] = 0.0f; + fftOut[i].real(0.0f); + fftOut[i].imag(0.0f); } } @@ -391,13 +452,13 @@ struct DataRxTx { } // calculate spectrum - std::copy(sampleAmplitudeAverage.begin(), sampleAmplitudeAverage.begin() + samplesPerFrame, fftIn); + std::copy(sampleAmplitudeAverage.begin(), sampleAmplitudeAverage.begin() + samplesPerFrame, fftIn.data()); - fftwf_execute(fftPlan); + FFT(fftIn.data(), fftOut.data(), samplesPerFrame, 1.0); double fsum = 0.0; for (int i = 0; i < samplesPerFrame; ++i) { - sampleSpectrum[i] = (fftOut[i][0]*fftOut[i][0] + fftOut[i][1]*fftOut[i][1]); + sampleSpectrum[i] = (fftOut[i].real()*fftOut[i].real() + fftOut[i].imag()*fftOut[i].imag()); fsum += sampleSpectrum[i]; } for (int i = 1; i < samplesPerFrame/2; ++i) { @@ -448,7 +509,7 @@ struct DataRxTx { std::copy( recordedAmplitude.begin() + offsetTx*step, - recordedAmplitude.begin() + offsetTx*step + samplesPerFrame, fftIn); + recordedAmplitude.begin() + offsetTx*step + samplesPerFrame, fftIn.data()); for (int k = 1; k < framesPerTx-1; ++k) { for (int i = 0; i < samplesPerFrame; ++i) { @@ -456,10 +517,10 @@ struct DataRxTx { } } - fftwf_execute(fftPlan); + FFT(fftIn.data(), fftOut.data(), samplesPerFrame, 1.0); for (int i = 0; i < samplesPerFrame; ++i) { - sampleSpectrum[i] = (fftOut[i][0]*fftOut[i][0] + fftOut[i][1]*fftOut[i][1]); + sampleSpectrum[i] = (fftOut[i].real()*fftOut[i].real() + fftOut[i].imag()*fftOut[i].imag()); } for (int i = 1; i < samplesPerFrame/2; ++i) { sampleSpectrum[i] += sampleSpectrum[samplesPerFrame - i]; @@ -638,9 +699,8 @@ struct DataRxTx { bool receivingData; bool analyzingData; - fftwf_plan fftPlan = 0; - float *fftIn; - fftwf_complex *fftOut = 0; + std::array fftIn; + std::array, kMaxSamplesPerFrame> fftOut; ::AmplitudeData sampleAmplitude; ::SpectrumData sampleSpectrum;