add beginning of DAB decoder and add missing hardware donors to the c…

…redits

CMakeLists.txt

@@ -43,6 +43,7 @@ option(OPT_BUILD_PORTAUDIO_SINK "Build PortAudio Sink Module (Dependencies: port
 
 # Decoders
 option(OPT_BUILD_ATV_DECODER "Build ATV decoder (no dependencies required)" OFF)
+option(OPT_BUILD_DAB_DECODER "Build the DAB/DAB+ decoder (no dependencies required)" OFF)
 option(OPT_BUILD_FALCON9_DECODER "Build the falcon9 live decoder (Dependencies: ffplay)" OFF)
 option(OPT_BUILD_KG_SSTV_DECODER "Build the KG SSTV (KG-STV) decoder module (no dependencies required)" OFF)
 option(OPT_BUILD_M17_DECODER "Build the M17 decoder module (Dependencies: codec2)" OFF)
@@ -242,6 +243,10 @@ if (OPT_BUILD_ATV_DECODER)
 add_subdirectory("decoder_modules/atv_decoder")
 endif (OPT_BUILD_ATV_DECODER)
 
+if (OPT_BUILD_DAB_DECODER)
+add_subdirectory("decoder_modules/dab_decoder")
+endif (OPT_BUILD_DAB_DECODER)
+
 if (OPT_BUILD_FALCON9_DECODER)
 add_subdirectory("decoder_modules/falcon9_decoder")
 endif (OPT_BUILD_FALCON9_DECODER)

core/src/credits.cpp

@@ -37,9 +37,12 @@ namespace sdrpp_credits {
     const char* hardwareDonators[] = {
         "Aaronia AG",
         "Airspy",
+        "Alex 4Z5LV",
         "Analog Devices",
         "CaribouLabs",
+        "Deepace",
         "Ettus Research",
+        "Harogic",
         "Howard Su",
         "MicroPhase",
         "Microtelecom",

decoder_modules/dab_decoder/CMakeLists.txt

@@ -0,0 +1,37 @@
+cmake_minimum_required(VERSION 3.13)
+project(dab_decoder)
+
+file(GLOB_RECURSE SRC "src/*.cpp" "src/*.c")
+
+include(${SDRPP_MODULE_CMAKE})
+
+target_include_directories(dab_decoder PRIVATE "src/")
+
+if (MSVC)
+    # Lib path
+    target_include_directories(dab_decoder PRIVATE "C:/Program Files/codec2/include/")
+    target_link_directories(dab_decoder PRIVATE "C:/Program Files/codec2/lib")
+
+    target_link_libraries(dab_decoder PRIVATE libcodec2)
+elseif (ANDROID)
+    target_include_directories(dab_decoder PUBLIC
+        /sdr-kit/${ANDROID_ABI}/include/codec2
+    )
+
+    target_link_libraries(dab_decoder PUBLIC
+        /sdr-kit/${ANDROID_ABI}/lib/libcodec2.so
+    )
+else ()
+    find_package(PkgConfig)
+
+    pkg_check_modules(LIBCODEC2 REQUIRED codec2)
+
+    target_include_directories(dab_decoder PRIVATE ${LIBCODEC2_INCLUDE_DIRS})
+    target_link_directories(dab_decoder PRIVATE ${LIBCODEC2_LIBRARY_DIRS})
+    target_link_libraries(dab_decoder PRIVATE ${LIBCODEC2_LIBRARIES})
+
+    # Include it because for some reason pkgconfig doesn't look here?
+    if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+        target_include_directories(dab_decoder PRIVATE "/usr/local/include")
+    endif()
+endif ()

decoder_modules/dab_decoder/src/dab_dsp.h

@@ -0,0 +1,280 @@
+#pragma once
+#include <dsp/processor.h>
+#include <utils/flog.h>
+#include <fftw3.h>
+#include "dab_phase_sym.h"
+
+namespace dab {
+    class CyclicSync : public dsp::Processor<dsp::complex_t, dsp::complex_t> {
+        using base_type = dsp::Processor<dsp::complex_t, dsp::complex_t>;
+    public:
+        CyclicSync() {}
+
+        // TODO: The default AGC rate is probably way too fast, plot out the avgCorr to see how much it moves
+        CyclicSync(dsp::stream<dsp::complex_t>* in, double symbolLength, double cyclicPrefixLength, double samplerate, float agcRate = 1e-3) { init(in, symbolLength, cyclicPrefixLength, samplerate, agcRate); }
+
+        void init(dsp::stream<dsp::complex_t>* in, double symbolLength, double cyclicPrefixLength, double samplerate, float agcRate = 1e-3) {
+            // Computer the number of samples for the symbol and its cyclic prefix
+            symbolSamps = round(samplerate * symbolLength);
+            prefixSamps = round(samplerate * cyclicPrefixLength);
+
+            // Allocate and clear the delay buffer
+            delayBuf = dsp::buffer::alloc<dsp::complex_t>(STREAM_BUFFER_SIZE + 64000);
+            dsp::buffer::clear(delayBuf, symbolSamps);
+
+            // Allocate and clear the history buffer
+            histBuf = dsp::buffer::alloc<dsp::complex_t>(prefixSamps);
+            dsp::buffer::clear(histBuf, prefixSamps);
+
+            // Compute the delay input addresses
+            delayBufInput = &delayBuf[symbolSamps];
+
+            // Compute the correlation AGC configuration
+            this->agcRate = agcRate;
+            agcRateInv = 1.0f - agcRate;
+
+            base_type::init(in);
+        }
+
+        void reset() {
+            assert(base_type::_block_init);
+            std::lock_guard<std::recursive_mutex> lck(base_type::ctrlMtx);
+            base_type::tempStop();
+
+            base_type::tempStart();
+        }
+
+        int run() {
+            int count = base_type::_in->read();
+            if (count < 0) { return -1; }
+
+            // Copy the data into the normal delay buffer
+            memcpy(delayBufInput, base_type::_in->readBuf, count * sizeof(dsp::complex_t));
+
+            // Flush the input stream
+            base_type::_in->flush();
+
+            // Do cross-correlation
+            for (int i = 0; i < count; i++) {
+                // Get the current history slot
+                dsp::complex_t* slot = &histBuf[histId++];
+
+                // Wrap around the history slot index (TODO: Check that the history buffer's length is correct)
+                histId %= prefixSamps;
+
+                // Kick out last value from the correlation
+                corr -= *slot;
+
+                // Save input value and compute the new prodct
+                dsp::complex_t val = delayBuf[i];
+                dsp::complex_t prod = val.conj()*delayBuf[i+symbolSamps];
+
+                // Add the new value to the correlation
+                *slot = prod;
+
+                // Add the new value to the history buffer
+                corr += prod;
+
+                // Compute sample amplitude
+                float rcorr = corr.amplitude();
+
+                // If a high enough peak is reached, reset the symbol counter
+                if (rcorr > avgCorr && rcorr > peakCorr) { // Note keeping an average level might not be needed
+                    peakCorr = rcorr;
+                    peakLCorr = lastCorr;
+                    samplesSincePeak = 0;
+                }
+
+                // If this is the sample right after the peak, save it
+                if (samplesSincePeak == 1) {
+                    peakRCorr = rcorr;
+                }
+
+                // Write the sample to the output
+                out.writeBuf[samplesSincePeak++] = val;
+
+                // If the end of the symbol is reached, send it off
+                if (samplesSincePeak >= symbolSamps) {
+                    if (!out.swap(symbolSamps)) {
+                        return -1;
+                    }
+                    samplesSincePeak = 0;
+                    peakCorr = 0;
+                }
+
+                // Update the average correlation
+                lastCorr = rcorr;
+
+                // Update the average correlation value
+                avgCorr = agcRate*rcorr + agcRateInv*avgCorr;
+            }
+
+            // Move unused data
+            memmove(delayBuf, &delayBuf[count], symbolSamps * sizeof(dsp::complex_t));
+
+            return count;
+        }
+
+    protected:
+        int symbolSamps;
+        int prefixSamps;
+
+        int histId = 0;
+        dsp::complex_t* histBuf;
+
+        dsp::complex_t* delayBuf;
+        dsp::complex_t* delayBufInput;
+
+        dsp::complex_t corr = { 0.0f, 0.0f };
+
+        int samplesSincePeak = 0;
+        float lastCorr = 0.0f;
+        float peakCorr = 0.0f;
+        float peakLCorr = 0.0f;
+        float peakRCorr = 0.0f;
+
+        // Note only required for DAB
+        float avgCorr = 0.0f;
+        float agcRate;
+        float agcRateInv;
+    };
+
+    class FrameFreqSync : public dsp::Processor<dsp::complex_t, dsp::complex_t> {
+        using base_type = dsp::Processor<dsp::complex_t, dsp::complex_t>;
+    public:
+        FrameFreqSync() {}
+
+        FrameFreqSync(dsp::stream<dsp::complex_t>* in, float agcRate = 0.01f) { init(in, agcRate); }
+
+        void init(dsp::stream<dsp::complex_t>* in, float agcRate = 0.01f) {
+            // Allocate buffers
+            amps = dsp::buffer::alloc<float>(2048);
+            conjRef = dsp::buffer::alloc<dsp::complex_t>(2048);
+            corrIn = (dsp::complex_t*)fftwf_alloc_complex(2048);
+            corrOut = (dsp::complex_t*)fftwf_alloc_complex(2048);
+
+            // Copy the phase reference
+            memcpy(conjRef, DAB_PHASE_SYM_CONJ, 2048 * sizeof(dsp::complex_t));
+
+            // Plan the FFT computation
+            plan = fftwf_plan_dft_1d(2048, (fftwf_complex*)corrIn, (fftwf_complex*)corrOut, FFTW_FORWARD, FFTW_ESTIMATE);
+
+            // Compute the correlation AGC configuration
+            this->agcRate = agcRate;
+            agcRateInv = 1.0f - agcRate;
+
+            base_type::init(in);
+        }
+
+        void reset() {
+            assert(base_type::_block_init);
+            std::lock_guard<std::recursive_mutex> lck(base_type::ctrlMtx);
+            base_type::tempStop();
+
+            base_type::tempStart();
+        }
+
+        int run() {
+            int count = base_type::_in->read();
+            if (count < 0) { return -1; }
+
+            // Apply frequency shift
+            lv_32fc_t phase = lv_cmake(1.0f, 0.0f);
+            lv_32fc_t phaseDelta = lv_cmake(cos(offset), sin(offset));
+#if VOLK_VERSION >= 030100
+            volk_32fc_s32fc_x2_rotator2_32fc((lv_32fc_t*)_in->readBuf, (lv_32fc_t*)_in->readBuf, phaseDelta, &phase, count);
+#else
+            volk_32fc_s32fc_x2_rotator_32fc((lv_32fc_t*)_in->readBuf, (lv_32fc_t*)_in->readBuf, phaseDelta, &phase, count);
+#endif
+
+            // Compute the amplitude amplitude of all samples
+            volk_32fc_magnitude_32f(amps, (lv_32fc_t*)_in->readBuf, 2048);
+
+            // Compute the average signal level by adding up all values
+            float level = 0.0f;
+            volk_32f_accumulator_s32f(&level, amps, 2048);
+
+            // Detect a frame sync condition
+            if (level < avgLvl * 0.5f) {
+                // Reset symbol counter
+                sym = 1;
+
+                // Update the average level
+                avgLvl = agcRate*level + agcRateInv*avgLvl;
+
+                // Flush the input stream and return
+                base_type::_in->flush();
+                return count;
+            }
+
+            // Update the average level
+            avgLvl = agcRate*level + agcRateInv*avgLvl;
+
+            // Handle phase reference
+            if (sym == 1) {
+                // Output the symbols (DEBUG ONLY)
+                memcpy(corrIn, _in->readBuf, 2048 * sizeof(dsp::complex_t));
+                fftwf_execute(plan);
+                volk_32fc_magnitude_32f(amps, (lv_32fc_t*)corrOut, 2048);
+                int outCount = 0;
+                dsp::complex_t pi4 = { cos(3.1415926535*0.25), sin(3.1415926535*0.25) };
+                for (int i = -767; i < 768; i++) {
+                    if (!i) { continue; }
+                    int cid0 = ((i-1) >= 0) ? (i-1) : 2048+(i-1);
+                    int cid1 = (i >= 0) ? i : 2048+i;;
+                    out.writeBuf[outCount++] = pi4 * (corrOut[cid1] * corrOut[cid0].conj()) * (1.0f/(amps[cid0]*amps[cid0]));
+                }
+                out.swap(outCount);
+
+                // Multiply the samples with the conjugated phase reference signal
+                volk_32fc_x2_multiply_32fc((lv_32fc_t*)corrIn, (lv_32fc_t*)_in->readBuf, (lv_32fc_t*)conjRef, 2048);
+
+                // Compute the FFT of the product
+                fftwf_execute(plan);
+
+                // Compute the amplitude of the bins
+                volk_32fc_magnitude_32f(amps, (lv_32fc_t*)corrOut, 2048);
+
+                // Locate highest power bin
+                uint32_t peakId;
+                volk_32f_index_max_32u(&peakId, amps, 2048);
+
+                // Obtain the value of the bins next to the peak
+                float peakL = amps[(peakId + 2047) % 2048];
+                float peakR = amps[(peakId + 1) % 2048];
+
+                // Compute the integer frequency offset
+                float offInt = (peakId < 1024) ? (float)peakId : ((float)peakId - 2048.0f);
+
+                // Compute the frequency offset in rad/samp
+                float off = 3.1415926535f * (offInt + ((peakR - peakL) / (peakR + peakL))) * (1.0f / 1024.0f);
+
+                // Run control loop
+                offset -= 0.1f*off;
+                flog::debug("Offset: {} Hz, Error: {} Hz, Avg Level: {}", offset * (0.5f/3.1415926535f)*2.048e6, off * (0.5f/3.1415926535f)*2.048e6, avgLvl);
+            }
+
+            // Increment the symbol counter
+            sym++;
+
+            // Flush the input stream and return
+            base_type::_in->flush();
+            return count;
+        }
+
+    protected:
+        fftwf_plan plan;
+
+        float* amps;
+        dsp::complex_t* conjRef;
+        dsp::complex_t* corrIn;
+        dsp::complex_t* corrOut;
+
+        int sym;
+        float offset = 0.0f;
+
+        float avgLvl = 0.0f;
+        float agcRate;
+        float agcRateInv;
+    };
+}