dsp(simd,resampler): optimize RTL widen, FIR scratch paths, and resampler layout

Author: arancormonk

include/dsd-neo/dsp/demod_state.h

@@ -83,8 +83,8 @@ struct demod_state {
     dsd_thread_t thread;
     float* lowpassed;
     double squelch_running_power;
-    float* resamp_taps; /* normalized taps, length = K*L */
-    float* resamp_hist; /* circular history, length = K */
+    float* resamp_taps; /* normalized taps as L contiguous phase blocks, length = K*L */
+    float* resamp_hist; /* mirrored history window, length = 2*K */
     int (*discriminator)(int, int, int, int);
     void (*mode_demod)(struct demod_state*);
     struct output_state* output_target;
@@ -167,7 +167,7 @@ struct demod_state {
     int resamp_phase;          /* 0..L-1 accumulator */
     int resamp_taps_len;       /* prototype taps length (padded to K*L) */
     int resamp_taps_per_phase; /* K = ceil(taps_len/L) */
-    int resamp_hist_head;      /* head index into circular history [0..K-1] */
+    int resamp_hist_head;      /* next write index into base history window [0..K-1] */
 
     /* Legacy FM FLL state (for non-CQPSK FM/C4FM paths).
      * Used by fll_update_error() and fll_mix_and_update() in demod_pipeline.cpp.

include/dsd-neo/dsp/resampler.h

@@ -23,9 +23,10 @@ struct demod_state;
 /**
  * @brief Design windowed-sinc low-pass prototype for polyphase upfirdn (runs at L*Fs_in).
  *
- * Taps are stored phase-major with stride L (k*L + phase). The function allocates
- * aligned storage for taps and history inside the provided demod_state and
- * initializes the resampler bookkeeping fields.
+ * Taps are stored as contiguous per-phase blocks with oldest-to-newest sample
+ * order inside each block. The function allocates aligned storage for taps and
+ * mirrored history inside the provided demod_state and initializes the
+ * resampler bookkeeping fields.
  *
  * @param s Demodulator state to receive resampler taps/history.
  * @param L Upsampling factor.

src/dsp/CMakeLists.txt

@@ -25,9 +25,11 @@ target_sources(dsd-neo_dsp PRIVATE
 # x86-64 SIMD sources with arch-specific flags
 if(CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|amd64")
   target_sources(dsd-neo_dsp PRIVATE simd_fir_sse2.cpp)
+  target_sources(dsd-neo_dsp PRIVATE simd_widen_sse2.cpp)
   # Explicit SSE2 flag for clarity and cross-compilation safety
   if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
     set_source_files_properties(simd_fir_sse2.cpp PROPERTIES COMPILE_FLAGS "-msse2")
+    set_source_files_properties(simd_widen_sse2.cpp PROPERTIES COMPILE_FLAGS "-msse2")
   endif()
 
   # AVX2+FMA needs explicit flag; add source only if the compiler accepts it
@@ -46,14 +48,17 @@ if(CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|amd64")
   endif()
   if(_avx2_flags)
     target_sources(dsd-neo_dsp PRIVATE simd_fir_avx2.cpp)
+    target_sources(dsd-neo_dsp PRIVATE simd_widen_avx2.cpp)
     set_source_files_properties(simd_fir_avx2.cpp PROPERTIES COMPILE_FLAGS "${_avx2_flags}")
+    set_source_files_properties(simd_widen_avx2.cpp PROPERTIES COMPILE_FLAGS "${_avx2_flags}")
     target_compile_definitions(dsd-neo_dsp PRIVATE $<$<COMPILE_LANGUAGE:CXX>:DSD_NEO_DSP_HAVE_AVX2_IMPL=1>)
   endif()
 endif()
 
 # ARM64 NEON (always available on AArch64)
 if(CMAKE_SYSTEM_PROCESSOR MATCHES "aarch64|arm64|ARM64")
   target_sources(dsd-neo_dsp PRIVATE simd_fir_neon.cpp)
+  target_sources(dsd-neo_dsp PRIVATE simd_widen_neon.cpp)
 endif()
 
 target_include_directories(dsd-neo_dsp

src/dsp/resampler.cpp

@@ -35,6 +35,7 @@
 #endif
 
 static const double kPi = 3.14159265358979323846;
+static const int kDefaultTapsPerPhase = 16;
 
 template <typename T>
 static inline T*
@@ -64,9 +65,37 @@ dsd_neo_sinc(double x) {
     return sin(kPi * x) / (kPi * x);
 }
 
+static inline float
+resamp_dot_contiguous(const float* DSD_NEO_RESTRICT samples, const float* DSD_NEO_RESTRICT taps, int len) {
+    float acc = 0.0f;
+    DSD_NEO_IVDEP
+    for (int k = 0; k < len; k++) {
+        acc += samples[k] * taps[k];
+    }
+    return acc;
+}
+
+static inline float
+resamp_dot16_contiguous(const float* DSD_NEO_RESTRICT samples, const float* DSD_NEO_RESTRICT taps) {
+    float acc0 = 0.0f;
+    float acc1 = 0.0f;
+    float acc2 = 0.0f;
+    float acc3 = 0.0f;
+
+    DSD_NEO_IVDEP
+    for (int k = 0; k < 16; k += 4) {
+        acc0 += samples[k + 0] * taps[k + 0];
+        acc1 += samples[k + 1] * taps[k + 1];
+        acc2 += samples[k + 2] * taps[k + 2];
+        acc3 += samples[k + 3] * taps[k + 3];
+    }
+
+    return (acc0 + acc1) + (acc2 + acc3);
+}
+
 void
 resamp_design(struct demod_state* s, int L, int M) {
-    int taps_per_phase = 16; /* K */
+    int taps_per_phase = kDefaultTapsPerPhase; /* K */
     if (taps_per_phase < 8) {
         taps_per_phase = 8;
     }
@@ -92,22 +121,22 @@ resamp_design(struct demod_state* s, int L, int M) {
         s->resamp_taps = (float*)mem_ptr;
     }
     {
-        void* mem_ptr = dsd_neo_aligned_malloc((size_t)taps_per_phase * sizeof(float));
+        void* mem_ptr = dsd_neo_aligned_malloc((size_t)taps_per_phase * 2U * sizeof(float));
         s->resamp_hist = (float*)mem_ptr;
     }
     if (!s->resamp_taps || !s->resamp_hist) {
         if (s->resamp_taps) {
-            free(s->resamp_taps);
+            dsd_neo_aligned_free(s->resamp_taps);
             s->resamp_taps = NULL;
         }
         if (s->resamp_hist) {
-            free(s->resamp_hist);
+            dsd_neo_aligned_free(s->resamp_hist);
             s->resamp_hist = NULL;
         }
         s->resamp_enabled = 0;
         return;
     }
-    memset(s->resamp_hist, 0, (size_t)taps_per_phase * sizeof(float));
+    memset(s->resamp_hist, 0, (size_t)taps_per_phase * 2U * sizeof(float));
     s->resamp_hist_head = 0;
 
     double gain = 0.0;
@@ -122,12 +151,16 @@ resamp_design(struct demod_state* s, int L, int M) {
     }
 
     const double phase_gain_comp = (double)L;
-    for (int n = 0; n < N; n++) {
-        int m = n - mid;
-        double w = 0.54 - 0.46 * cos(2.0 * kPi * (double)n / (double)(N - 1));
-        double h = 2.0 * fc * dsd_neo_sinc(2.0 * fc * (double)m);
-        double t = (h * w / gain) * phase_gain_comp;
-        s->resamp_taps[n] = (float)t;
+    for (int phase = 0; phase < L; phase++) {
+        float* phase_taps = s->resamp_taps + (size_t)phase * (size_t)taps_per_phase;
+        for (int k = 0; k < taps_per_phase; k++) {
+            int src_index = phase + ((taps_per_phase - 1 - k) * L);
+            int m = src_index - mid;
+            double w = 0.54 - 0.46 * cos(2.0 * kPi * (double)src_index / (double)(N - 1));
+            double h = 2.0 * fc * dsd_neo_sinc(2.0 * fc * (double)m);
+            double t = (h * w / gain) * phase_gain_comp;
+            phase_taps[k] = (float)t;
+        }
     }
 
     s->resamp_L = L;
@@ -159,26 +192,17 @@ resamp_process_block(struct demod_state* s, const float* DSD_NEO_RESTRICT in, in
 
     for (int n = 0; n < in_len; n++) {
         hist[head] = in_al[n];
+        hist[head + K] = in_al[n];
         head++;
         if (head == K) {
             head = 0;
         }
         int local_phase = phase;
+        const float* DSD_NEO_RESTRICT hist_window = hist + head;
         while (local_phase < L) {
-            float acc = 0.0f;
-            const float* DSD_NEO_RESTRICT tk = taps_al + local_phase;
-            int idx = head - 1;
-            if (idx < 0) {
-                idx += K;
-            }
-            for (int k = 0; k < K; k++) {
-                acc += hist[idx] * tk[0];
-                tk += L;
-                idx--;
-                if (idx < 0) {
-                    idx += K;
-                }
-            }
+            const float* DSD_NEO_RESTRICT phase_taps = taps_al + (size_t)local_phase * (size_t)K;
+            float acc = (K == kDefaultTapsPerPhase) ? resamp_dot16_contiguous(hist_window, phase_taps)
+                                                    : resamp_dot_contiguous(hist_window, phase_taps, K);
             out_al[out_len++] = acc;
             local_phase += M;
         }

src/dsp/simd_fir.cpp

@@ -16,24 +16,7 @@
 #include <atomic>
 #include <cstring>
 
-/* Platform-specific CPU feature detection for optional AVX2 dispatch */
-#if defined(__x86_64__) || defined(_M_X64)
-#if defined(DSD_NEO_DSP_HAVE_AVX2_IMPL)
-#if defined(_MSC_VER)
-#include <intrin.h>
-#else
-#include <cpuid.h>
-
-/* Use inline assembly for _xgetbv to avoid target-specific option issues */
-static inline unsigned long long
-dsd_xgetbv(unsigned int xcr) {
-    unsigned int eax, edx;
-    __asm__ __volatile__("xgetbv" : "=a"(eax), "=d"(edx) : "c"(xcr));
-    return ((unsigned long long)edx << 32) | eax;
-}
-#endif
-#endif
-#endif
+#include "simd_x86_cpu.h"
 
 /* Forward declarations for SIMD specializations (defined in arch-specific TUs) */
 #if defined(__x86_64__) || defined(_M_X64)
@@ -305,61 +288,6 @@ simd_hb_decim2_real_scalar(const float* in, int in_len, float* out, float* hist,
     return out_len;
 }
 
-/* -------------------------------------------------------------------------- */
-/* CPU Feature Detection                                                      */
-/* -------------------------------------------------------------------------- */
-
-#if (defined(__x86_64__) || defined(_M_X64)) && defined(DSD_NEO_DSP_HAVE_AVX2_IMPL)
-
-static bool
-cpu_has_avx2_with_os_support() {
-#if defined(__GNUC__) || defined(__clang__)
-    unsigned int eax, ebx, ecx, edx;
-    if (!__get_cpuid(1, &eax, &ebx, &ecx, &edx)) {
-        return false;
-    }
-    bool osxsave = (ecx & (1u << 27)) != 0;
-    bool avx = (ecx & (1u << 28)) != 0;
-    bool fma = (ecx & (1u << 12)) != 0;
-    if (!osxsave || !avx) {
-        return false;
-    }
-    /* Check OS has enabled YMM state saving via XGETBV */
-    unsigned long long xcr0 = dsd_xgetbv(0);
-    bool ymm_enabled = (xcr0 & 0x6) == 0x6; /* XMM + YMM state enabled */
-    if (!ymm_enabled) {
-        return false;
-    }
-    /* Check AVX2 support in extended features */
-    if (!__get_cpuid_count(7, 0, &eax, &ebx, &ecx, &edx)) {
-        return false;
-    }
-    bool has_avx2 = (ebx & (1u << 5)) != 0;
-    return has_avx2 && fma;
-#elif defined(_MSC_VER)
-    int cpuInfo[4];
-    __cpuid(cpuInfo, 1);
-    bool osxsave = (cpuInfo[2] & (1 << 27)) != 0;
-    bool avx = (cpuInfo[2] & (1 << 28)) != 0;
-    bool fma = (cpuInfo[2] & (1 << 12)) != 0;
-    if (!osxsave || !avx) {
-        return false;
-    }
-    unsigned long long xcr0 = _xgetbv(0);
-    bool ymm_enabled = (xcr0 & 0x6) == 0x6;
-    if (!ymm_enabled) {
-        return false;
-    }
-    __cpuidex(cpuInfo, 7, 0);
-    bool has_avx2 = (cpuInfo[1] & (1 << 5)) != 0;
-    return has_avx2 && fma;
-#else
-    return false;
-#endif
-}
-
-#endif /* x86_64 && DSD_NEO_DSP_HAVE_AVX2_IMPL */
-
 /* -------------------------------------------------------------------------- */
 /* Function Pointer Dispatch                                                  */
 /* -------------------------------------------------------------------------- */
@@ -390,7 +318,7 @@ simd_fir_init_dispatch() {
     /* Perform one-time initialization */
 #if defined(__x86_64__) || defined(_M_X64)
 #if defined(DSD_NEO_DSP_HAVE_AVX2_IMPL)
-    if (cpu_has_avx2_with_os_support()) {
+    if (dsd_neo_cpu_has_avx2_with_os_support()) {
         g_fir_complex_impl = simd_fir_complex_apply_avx2;
         g_hb_decim2_complex_impl = simd_hb_decim2_complex_avx2;
         g_hb_decim2_real_impl = simd_hb_decim2_real_avx2;