DSPark 1.6.1
Header-only audio DSP framework in pure C++20 — zero dependencies
Loading...
Searching...
No Matches
SimdOps.h
Go to the documentation of this file.
1// DSPark -- Professional Audio DSP Framework
2// Copyright (c) 2026 Cristian Moresi -- MIT License
3
4#pragma once
5
20// --- Platform SIMD detection ------------------------------------------------
21#if defined(_M_AMD64) || defined(_M_X64) || defined(__x86_64__) || defined(__amd64__)
22 #define DSPARK_SIMD_SSE2 1
23 #include <emmintrin.h> // SSE2
24 #if defined(__AVX__)
25 #define DSPARK_SIMD_AVX 1
26 #include <immintrin.h> // AVX
27 #endif
28 // Detect FMA3 support. __FMA__ is the authoritative flag on GCC/Clang:
29 // AVX2 alone does not licence FMA codegen there, and calling _mm*_fmadd_*
30 // without the fma target feature is a hard compile error (e.g. building
31 // with -mavx2 but not -mfma). MSVC never defines __FMA__ but permits FMA
32 // intrinsics under /arch:AVX2; clang-cl defines _MSC_VER yet enforces
33 // target features like clang, so it must take the __FMA__ route only.
34 #if defined(__FMA__) || (defined(__AVX2__) && defined(_MSC_VER) && !defined(__clang__))
35 #define DSPARK_SIMD_FMA 1
36 #endif
37#elif defined(__aarch64__) || defined(_M_ARM64)
38 #define DSPARK_SIMD_NEON 1
39 #include <arm_neon.h>
40#elif defined(__EMSCRIPTEN__) && defined(__wasm_simd128__) && defined(__SSE2__)
41 // Emscripten lowers the SSE2 intrinsic set to Wasm SIMD128 when compiled
42 // with -msimd128 -msse2 (both flags required). Anything less stays scalar.
43 #define DSPARK_SIMD_SSE2 1
44 #include <emmintrin.h>
45#endif
46
47// --- Restrict Macro for Pointer Aliasing Optimization -----------------------
48// NOTE: restrict constrains WRITTEN ranges only. A range written through one
49// DSPARK_RESTRICT pointer must not overlap any other range accessed in the
50// same call (so dst != src, no partial overlaps). Read-only pointers may
51// alias each other freely: dotProduct(data, data, n) is well defined, and is
52// exactly what sumOfSquares() does. Defined idempotently (#ifndef) because a
53// few headers historically declared the same macro.
54#ifndef DSPARK_RESTRICT
55 #if defined(_MSC_VER)
56 #define DSPARK_RESTRICT __restrict
57 #elif defined(__GNUC__) || defined(__clang__)
58 #define DSPARK_RESTRICT __restrict__
59 #else
60 #define DSPARK_RESTRICT
61 #endif
62#endif
63
64#include <cstdint>
65#include <type_traits>
66
67namespace dspark {
68namespace simd {
69
70// ============================================================================
71// addWithGain -- dst[i] += src[i] * gain
72// ============================================================================
73
77inline void addWithGain(float* DSPARK_RESTRICT dst, const float* DSPARK_RESTRICT src, float gain, int count) noexcept
78{
79#if defined(DSPARK_SIMD_AVX)
80 const __m256 vGain = _mm256_set1_ps(gain);
81 int i = 0;
82 for (; i + 7 < count; i += 8)
83 {
84 __m256 vDst = _mm256_loadu_ps(dst + i);
85 __m256 vSrc = _mm256_loadu_ps(src + i);
86
87 #if defined(DSPARK_SIMD_FMA)
88 _mm256_storeu_ps(dst + i, _mm256_fmadd_ps(vSrc, vGain, vDst));
89 #else
90 _mm256_storeu_ps(dst + i, _mm256_add_ps(vDst, _mm256_mul_ps(vSrc, vGain)));
91 #endif
92 }
93 for (; i < count; ++i) dst[i] += src[i] * gain;
94
95#elif defined(DSPARK_SIMD_SSE2)
96 const __m128 vGain = _mm_set1_ps(gain);
97 int i = 0;
98 for (; i + 3 < count; i += 4)
99 {
100 __m128 vDst = _mm_loadu_ps(dst + i);
101 __m128 vSrc = _mm_loadu_ps(src + i);
102
103 #if defined(DSPARK_SIMD_FMA)
104 _mm_storeu_ps(dst + i, _mm_fmadd_ps(vSrc, vGain, vDst));
105 #else
106 _mm_storeu_ps(dst + i, _mm_add_ps(vDst, _mm_mul_ps(vSrc, vGain)));
107 #endif
108 }
109 for (; i < count; ++i) dst[i] += src[i] * gain;
110
111#elif defined(DSPARK_SIMD_NEON)
112 const float32x4_t vGain = vdupq_n_f32(gain);
113 int i = 0;
114 for (; i + 3 < count; i += 4)
115 {
116 float32x4_t vDst = vld1q_f32(dst + i);
117 float32x4_t vSrc = vld1q_f32(src + i);
118 // vfmaq guarantees a fused FMLA; vmlaq has separate mul+add semantics
119 // in ACLE and GCC lowers it to fmul+fadd.
120 vst1q_f32(dst + i, vfmaq_f32(vDst, vSrc, vGain));
121 }
122 for (; i < count; ++i) dst[i] += src[i] * gain;
123
124#else
125 for (int i = 0; i < count; ++i) dst[i] += src[i] * gain;
126#endif
127}
128
130inline void addWithGain(double* DSPARK_RESTRICT dst, const double* DSPARK_RESTRICT src, double gain, int count) noexcept
131{
132#if defined(DSPARK_SIMD_AVX)
133 const __m256d vGain = _mm256_set1_pd(gain);
134 int i = 0;
135 for (; i + 3 < count; i += 4)
136 {
137 __m256d vDst = _mm256_loadu_pd(dst + i);
138 __m256d vSrc = _mm256_loadu_pd(src + i);
139
140 #if defined(DSPARK_SIMD_FMA)
141 _mm256_storeu_pd(dst + i, _mm256_fmadd_pd(vSrc, vGain, vDst));
142 #else
143 _mm256_storeu_pd(dst + i, _mm256_add_pd(vDst, _mm256_mul_pd(vSrc, vGain)));
144 #endif
145 }
146 for (; i < count; ++i) dst[i] += src[i] * gain;
147
148#elif defined(DSPARK_SIMD_SSE2)
149 const __m128d vGain = _mm_set1_pd(gain);
150 int i = 0;
151 for (; i + 1 < count; i += 2)
152 {
153 __m128d vDst = _mm_loadu_pd(dst + i);
154 __m128d vSrc = _mm_loadu_pd(src + i);
155
156 #if defined(DSPARK_SIMD_FMA)
157 _mm_storeu_pd(dst + i, _mm_fmadd_pd(vSrc, vGain, vDst));
158 #else
159 _mm_storeu_pd(dst + i, _mm_add_pd(vDst, _mm_mul_pd(vSrc, vGain)));
160 #endif
161 }
162 for (; i < count; ++i) dst[i] += src[i] * gain;
163
164#elif defined(DSPARK_SIMD_NEON)
165 const float64x2_t vGain = vdupq_n_f64(gain);
166 int i = 0;
167 for (; i + 1 < count; i += 2)
168 {
169 float64x2_t vDst = vld1q_f64(dst + i);
170 float64x2_t vSrc = vld1q_f64(src + i);
171 vst1q_f64(dst + i, vfmaq_f64(vDst, vSrc, vGain));
172 }
173 for (; i < count; ++i) dst[i] += src[i] * gain;
174#else
175 for (int i = 0; i < count; ++i) dst[i] += src[i] * gain;
176#endif
177}
178
179// ============================================================================
180// applyGain -- data[i] *= gain
181// ============================================================================
182
186inline void applyGain(float* DSPARK_RESTRICT data, float gain, int count) noexcept
187{
188#if defined(DSPARK_SIMD_AVX)
189 const __m256 vGain = _mm256_set1_ps(gain);
190 int i = 0;
191 for (; i + 7 < count; i += 8)
192 {
193 __m256 v = _mm256_loadu_ps(data + i);
194 _mm256_storeu_ps(data + i, _mm256_mul_ps(v, vGain));
195 }
196 for (; i < count; ++i) data[i] *= gain;
197
198#elif defined(DSPARK_SIMD_SSE2)
199 const __m128 vGain = _mm_set1_ps(gain);
200 int i = 0;
201 for (; i + 3 < count; i += 4)
202 {
203 __m128 v = _mm_loadu_ps(data + i);
204 _mm_storeu_ps(data + i, _mm_mul_ps(v, vGain));
205 }
206 for (; i < count; ++i) data[i] *= gain;
207
208#elif defined(DSPARK_SIMD_NEON)
209 const float32x4_t vGain = vdupq_n_f32(gain);
210 int i = 0;
211 for (; i + 3 < count; i += 4)
212 {
213 float32x4_t v = vld1q_f32(data + i);
214 vst1q_f32(data + i, vmulq_f32(v, vGain));
215 }
216 for (; i < count; ++i) data[i] *= gain;
217
218#else
219 for (int i = 0; i < count; ++i) data[i] *= gain;
220#endif
221}
222
224inline void applyGain(double* DSPARK_RESTRICT data, double gain, int count) noexcept
225{
226#if defined(DSPARK_SIMD_AVX)
227 const __m256d vGain = _mm256_set1_pd(gain);
228 int i = 0;
229 for (; i + 3 < count; i += 4)
230 _mm256_storeu_pd(data + i, _mm256_mul_pd(_mm256_loadu_pd(data + i), vGain));
231 for (; i < count; ++i) data[i] *= gain;
232
233#elif defined(DSPARK_SIMD_SSE2)
234 const __m128d vGain = _mm_set1_pd(gain);
235 int i = 0;
236 for (; i + 1 < count; i += 2)
237 {
238 __m128d v = _mm_loadu_pd(data + i);
239 _mm_storeu_pd(data + i, _mm_mul_pd(v, vGain));
240 }
241 for (; i < count; ++i) data[i] *= gain;
242
243#elif defined(DSPARK_SIMD_NEON)
244 const float64x2_t vGain = vdupq_n_f64(gain);
245 int i = 0;
246 for (; i + 1 < count; i += 2)
247 vst1q_f64(data + i, vmulq_f64(vld1q_f64(data + i), vGain));
248 for (; i < count; ++i) data[i] *= gain;
249#else
250 for (int i = 0; i < count; ++i) data[i] *= gain;
251#endif
252}
253
254// ============================================================================
255// peakLevel -- max(abs(data[i]))
256// ============================================================================
257
264inline float peakLevel(const float* DSPARK_RESTRICT data, int count) noexcept
265{
266#if defined(DSPARK_SIMD_AVX)
267 const __m256 absMask = _mm256_castsi256_ps(_mm256_set1_epi32(0x7FFFFFFF));
268 __m256 vMax = _mm256_setzero_ps();
269 int i = 0;
270 for (; i + 7 < count; i += 8)
271 {
272 __m256 v = _mm256_loadu_ps(data + i);
273 v = _mm256_and_ps(v, absMask);
274 // Data operand first: MAXPS returns the SECOND operand on NaN, so a
275 // NaN sample is skipped and the accumulator keeps its history
276 // (matching the scalar tail, whose comparison also skips NaNs).
277 vMax = _mm256_max_ps(v, vMax);
278 }
279 __m128 hi = _mm256_extractf128_ps(vMax, 1);
280 __m128 lo = _mm256_castps256_ps128(vMax);
281 __m128 m4 = _mm_max_ps(lo, hi);
282 __m128 m2 = _mm_max_ps(m4, _mm_movehl_ps(m4, m4));
283 __m128 m1 = _mm_max_ss(m2, _mm_shuffle_ps(m2, m2, 1));
284 float peak = _mm_cvtss_f32(m1);
285
286 for (; i < count; ++i)
287 {
288 float a = data[i] < 0.0f ? -data[i] : data[i];
289 if (a > peak) peak = a;
290 }
291 return peak;
292
293#elif defined(DSPARK_SIMD_SSE2)
294 const __m128 absMask = _mm_castsi128_ps(_mm_set1_epi32(0x7FFFFFFF));
295 __m128 vMax = _mm_setzero_ps();
296 int i = 0;
297 for (; i + 3 < count; i += 4)
298 {
299 __m128 v = _mm_loadu_ps(data + i);
300 v = _mm_and_ps(v, absMask);
301 vMax = _mm_max_ps(v, vMax); // data first: NaN skipped (see AVX note)
302 }
303 __m128 shuf = _mm_movehl_ps(vMax, vMax);
304 __m128 maxPair = _mm_max_ps(vMax, shuf);
305 __m128 maxSingle = _mm_max_ss(maxPair, _mm_shuffle_ps(maxPair, maxPair, 1));
306 float peak = _mm_cvtss_f32(maxSingle);
307
308 for (; i < count; ++i)
309 {
310 float a = data[i] < 0.0f ? -data[i] : data[i];
311 if (a > peak) peak = a;
312 }
313 return peak;
314
315#elif defined(DSPARK_SIMD_NEON)
316 float32x4_t vMax = vdupq_n_f32(0.0f);
317 int i = 0;
318 for (; i + 3 < count; i += 4)
319 {
320 float32x4_t v = vld1q_f32(data + i);
321 v = vabsq_f32(v);
322 vMax = vmaxnmq_f32(vMax, v); // FMAXNM: NaN operands are skipped
323 }
324 float peak = vmaxvq_f32(vMax);
325
326 for (; i < count; ++i)
327 {
328 float a = data[i] < 0.0f ? -data[i] : data[i];
329 if (a > peak) peak = a;
330 }
331 return peak;
332
333#else
334 float peak = 0.0f;
335 for (int i = 0; i < count; ++i)
336 {
337 float a = data[i] < 0.0f ? -data[i] : data[i];
338 if (a > peak) peak = a;
339 }
340 return peak;
341#endif
342}
343
345inline double peakLevel(const double* DSPARK_RESTRICT data, int count) noexcept
346{
347#if defined(DSPARK_SIMD_AVX)
348 const __m256d absMask = _mm256_castsi256_pd(_mm256_set1_epi64x(0x7FFFFFFFFFFFFFFFll));
349 __m256d vMax = _mm256_setzero_pd();
350 int i = 0;
351 for (; i + 3 < count; i += 4)
352 vMax = _mm256_max_pd(_mm256_and_pd(_mm256_loadu_pd(data + i), absMask), vMax);
353 __m128d lo = _mm256_castpd256_pd128(vMax);
354 __m128d hi = _mm256_extractf128_pd(vMax, 1);
355 __m128d m2 = _mm_max_pd(lo, hi);
356 double peak = _mm_cvtsd_f64(_mm_max_sd(m2, _mm_unpackhi_pd(m2, m2)));
357
358 for (; i < count; ++i)
359 {
360 double a = data[i] < 0.0 ? -data[i] : data[i];
361 if (a > peak) peak = a;
362 }
363 return peak;
364
365#elif defined(DSPARK_SIMD_SSE2)
366 const __m128d absMask = _mm_castsi128_pd(_mm_set_epi64x(
367 static_cast<int64_t>(0x7FFFFFFFFFFFFFFF),
368 static_cast<int64_t>(0x7FFFFFFFFFFFFFFF)));
369 __m128d vMax = _mm_setzero_pd();
370 int i = 0;
371 for (; i + 1 < count; i += 2)
372 {
373 __m128d v = _mm_loadu_pd(data + i);
374 v = _mm_and_pd(v, absMask);
375 vMax = _mm_max_pd(v, vMax); // data first: NaN skipped
376 }
377 __m128d hi = _mm_unpackhi_pd(vMax, vMax);
378 __m128d maxVal = _mm_max_sd(vMax, hi);
379 double peak = _mm_cvtsd_f64(maxVal);
380
381 for (; i < count; ++i)
382 {
383 double a = data[i] < 0.0 ? -data[i] : data[i];
384 if (a > peak) peak = a;
385 }
386 return peak;
387
388#elif defined(DSPARK_SIMD_NEON)
389 float64x2_t vMax = vdupq_n_f64(0.0);
390 int i = 0;
391 for (; i + 1 < count; i += 2)
392 vMax = vmaxnmq_f64(vMax, vabsq_f64(vld1q_f64(data + i)));
393 double peak = vmaxvq_f64(vMax);
394
395 for (; i < count; ++i)
396 {
397 double a = data[i] < 0.0 ? -data[i] : data[i];
398 if (a > peak) peak = a;
399 }
400 return peak;
401#else
402 double peak = 0.0;
403 for (int i = 0; i < count; ++i)
404 {
405 double a = data[i] < 0.0 ? -data[i] : data[i];
406 if (a > peak) peak = a;
407 }
408 return peak;
409#endif
410}
411
412// ============================================================================
413// dotProduct -- sum(a[i] * b[i])
414// ============================================================================
415
419inline float dotProduct(const float* DSPARK_RESTRICT a, const float* DSPARK_RESTRICT b, int count) noexcept
420{
421#if defined(DSPARK_SIMD_AVX)
422 // Four independent accumulators break the FMA latency chain (4-5 cycles),
423 // letting long FIR kernels run at full multiply throughput (~2-4x faster).
424 __m256 vSum0 = _mm256_setzero_ps();
425 __m256 vSum1 = _mm256_setzero_ps();
426 __m256 vSum2 = _mm256_setzero_ps();
427 __m256 vSum3 = _mm256_setzero_ps();
428 int i = 0;
429 for (; i + 31 < count; i += 32)
430 {
431 #if defined(DSPARK_SIMD_FMA)
432 vSum0 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i), _mm256_loadu_ps(b + i), vSum0);
433 vSum1 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i + 8), _mm256_loadu_ps(b + i + 8), vSum1);
434 vSum2 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i + 16), _mm256_loadu_ps(b + i + 16), vSum2);
435 vSum3 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i + 24), _mm256_loadu_ps(b + i + 24), vSum3);
436 #else
437 vSum0 = _mm256_add_ps(vSum0, _mm256_mul_ps(_mm256_loadu_ps(a + i), _mm256_loadu_ps(b + i)));
438 vSum1 = _mm256_add_ps(vSum1, _mm256_mul_ps(_mm256_loadu_ps(a + i + 8), _mm256_loadu_ps(b + i + 8)));
439 vSum2 = _mm256_add_ps(vSum2, _mm256_mul_ps(_mm256_loadu_ps(a + i + 16), _mm256_loadu_ps(b + i + 16)));
440 vSum3 = _mm256_add_ps(vSum3, _mm256_mul_ps(_mm256_loadu_ps(a + i + 24), _mm256_loadu_ps(b + i + 24)));
441 #endif
442 }
443 for (; i + 7 < count; i += 8)
444 {
445 #if defined(DSPARK_SIMD_FMA)
446 vSum0 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i), _mm256_loadu_ps(b + i), vSum0);
447 #else
448 vSum0 = _mm256_add_ps(vSum0, _mm256_mul_ps(_mm256_loadu_ps(a + i), _mm256_loadu_ps(b + i)));
449 #endif
450 }
451 __m256 vSum = _mm256_add_ps(_mm256_add_ps(vSum0, vSum1), _mm256_add_ps(vSum2, vSum3));
452 __m128 hi = _mm256_extractf128_ps(vSum, 1);
453 __m128 lo = _mm256_castps256_ps128(vSum);
454 __m128 sum4 = _mm_add_ps(lo, hi);
455 __m128 shuf = _mm_movehl_ps(sum4, sum4);
456 __m128 sum2 = _mm_add_ps(sum4, shuf);
457 __m128 sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
458 float sum = _mm_cvtss_f32(sum1);
459
460 for (; i < count; ++i) sum += a[i] * b[i];
461 return sum;
462
463#elif defined(DSPARK_SIMD_SSE2)
464 // Two accumulators hide the add/FMA latency on SSE-class hardware.
465 __m128 vSum0 = _mm_setzero_ps();
466 __m128 vSum1 = _mm_setzero_ps();
467 int i = 0;
468 for (; i + 7 < count; i += 8)
469 {
470 #if defined(DSPARK_SIMD_FMA)
471 vSum0 = _mm_fmadd_ps(_mm_loadu_ps(a + i), _mm_loadu_ps(b + i), vSum0);
472 vSum1 = _mm_fmadd_ps(_mm_loadu_ps(a + i + 4), _mm_loadu_ps(b + i + 4), vSum1);
473 #else
474 vSum0 = _mm_add_ps(vSum0, _mm_mul_ps(_mm_loadu_ps(a + i), _mm_loadu_ps(b + i)));
475 vSum1 = _mm_add_ps(vSum1, _mm_mul_ps(_mm_loadu_ps(a + i + 4), _mm_loadu_ps(b + i + 4)));
476 #endif
477 }
478 for (; i + 3 < count; i += 4)
479 {
480 #if defined(DSPARK_SIMD_FMA)
481 vSum0 = _mm_fmadd_ps(_mm_loadu_ps(a + i), _mm_loadu_ps(b + i), vSum0);
482 #else
483 vSum0 = _mm_add_ps(vSum0, _mm_mul_ps(_mm_loadu_ps(a + i), _mm_loadu_ps(b + i)));
484 #endif
485 }
486 alignas(16) float tmp[4];
487 _mm_store_ps(tmp, _mm_add_ps(vSum0, vSum1));
488 float sum = tmp[0] + tmp[1] + tmp[2] + tmp[3];
489
490 for (; i < count; ++i) sum += a[i] * b[i];
491 return sum;
492
493#elif defined(DSPARK_SIMD_NEON)
494 float32x4_t vSum0 = vdupq_n_f32(0.0f);
495 float32x4_t vSum1 = vdupq_n_f32(0.0f);
496 int i = 0;
497 for (; i + 7 < count; i += 8)
498 {
499 vSum0 = vfmaq_f32(vSum0, vld1q_f32(a + i), vld1q_f32(b + i));
500 vSum1 = vfmaq_f32(vSum1, vld1q_f32(a + i + 4), vld1q_f32(b + i + 4));
501 }
502 for (; i + 3 < count; i += 4)
503 vSum0 = vfmaq_f32(vSum0, vld1q_f32(a + i), vld1q_f32(b + i));
504 float sum = vaddvq_f32(vaddq_f32(vSum0, vSum1));
505
506 for (; i < count; ++i) sum += a[i] * b[i];
507 return sum;
508
509#else
510 float sum = 0.0f;
511 for (int i = 0; i < count; ++i) sum += a[i] * b[i];
512 return sum;
513#endif
514}
515
517inline double dotProduct(const double* DSPARK_RESTRICT a, const double* DSPARK_RESTRICT b, int count) noexcept
518{
519#if defined(DSPARK_SIMD_AVX)
520 __m256d vSum0 = _mm256_setzero_pd();
521 __m256d vSum1 = _mm256_setzero_pd();
522 int i = 0;
523 for (; i + 7 < count; i += 8)
524 {
525 #if defined(DSPARK_SIMD_FMA)
526 vSum0 = _mm256_fmadd_pd(_mm256_loadu_pd(a + i), _mm256_loadu_pd(b + i), vSum0);
527 vSum1 = _mm256_fmadd_pd(_mm256_loadu_pd(a + i + 4), _mm256_loadu_pd(b + i + 4), vSum1);
528 #else
529 vSum0 = _mm256_add_pd(vSum0, _mm256_mul_pd(_mm256_loadu_pd(a + i), _mm256_loadu_pd(b + i)));
530 vSum1 = _mm256_add_pd(vSum1, _mm256_mul_pd(_mm256_loadu_pd(a + i + 4), _mm256_loadu_pd(b + i + 4)));
531 #endif
532 }
533 for (; i + 3 < count; i += 4)
534 {
535 #if defined(DSPARK_SIMD_FMA)
536 vSum0 = _mm256_fmadd_pd(_mm256_loadu_pd(a + i), _mm256_loadu_pd(b + i), vSum0);
537 #else
538 vSum0 = _mm256_add_pd(vSum0, _mm256_mul_pd(_mm256_loadu_pd(a + i), _mm256_loadu_pd(b + i)));
539 #endif
540 }
541 __m256d vSum = _mm256_add_pd(vSum0, vSum1);
542 __m128d lo = _mm256_castpd256_pd128(vSum);
543 __m128d hi = _mm256_extractf128_pd(vSum, 1);
544 __m128d s2 = _mm_add_pd(lo, hi);
545 double sum = _mm_cvtsd_f64(_mm_add_sd(s2, _mm_unpackhi_pd(s2, s2)));
546
547 for (; i < count; ++i) sum += a[i] * b[i];
548 return sum;
549
550#elif defined(DSPARK_SIMD_SSE2)
551 __m128d vSum0 = _mm_setzero_pd();
552 __m128d vSum1 = _mm_setzero_pd();
553 int i = 0;
554 for (; i + 3 < count; i += 4)
555 {
556 #if defined(DSPARK_SIMD_FMA)
557 vSum0 = _mm_fmadd_pd(_mm_loadu_pd(a + i), _mm_loadu_pd(b + i), vSum0);
558 vSum1 = _mm_fmadd_pd(_mm_loadu_pd(a + i + 2), _mm_loadu_pd(b + i + 2), vSum1);
559 #else
560 vSum0 = _mm_add_pd(vSum0, _mm_mul_pd(_mm_loadu_pd(a + i), _mm_loadu_pd(b + i)));
561 vSum1 = _mm_add_pd(vSum1, _mm_mul_pd(_mm_loadu_pd(a + i + 2), _mm_loadu_pd(b + i + 2)));
562 #endif
563 }
564 for (; i + 1 < count; i += 2)
565 {
566 #if defined(DSPARK_SIMD_FMA)
567 vSum0 = _mm_fmadd_pd(_mm_loadu_pd(a + i), _mm_loadu_pd(b + i), vSum0);
568 #else
569 vSum0 = _mm_add_pd(vSum0, _mm_mul_pd(_mm_loadu_pd(a + i), _mm_loadu_pd(b + i)));
570 #endif
571 }
572 __m128d vSum = _mm_add_pd(vSum0, vSum1);
573 __m128d hi = _mm_unpackhi_pd(vSum, vSum);
574 double sum = _mm_cvtsd_f64(_mm_add_sd(vSum, hi));
575
576 for (; i < count; ++i) sum += a[i] * b[i];
577 return sum;
578
579#elif defined(DSPARK_SIMD_NEON)
580 float64x2_t vSum0 = vdupq_n_f64(0.0);
581 float64x2_t vSum1 = vdupq_n_f64(0.0);
582 int i = 0;
583 for (; i + 3 < count; i += 4)
584 {
585 vSum0 = vfmaq_f64(vSum0, vld1q_f64(a + i), vld1q_f64(b + i));
586 vSum1 = vfmaq_f64(vSum1, vld1q_f64(a + i + 2), vld1q_f64(b + i + 2));
587 }
588 for (; i + 1 < count; i += 2)
589 vSum0 = vfmaq_f64(vSum0, vld1q_f64(a + i), vld1q_f64(b + i));
590 double sum = vaddvq_f64(vaddq_f64(vSum0, vSum1));
591
592 for (; i < count; ++i) sum += a[i] * b[i];
593 return sum;
594#else
595 double sum = 0.0;
596 for (int i = 0; i < count; ++i) sum += a[i] * b[i];
597 return sum;
598#endif
599}
600
601// ============================================================================
602// add -- dst[i] += src[i]
603// ============================================================================
604
608inline void add(float* DSPARK_RESTRICT dst, const float* DSPARK_RESTRICT src, int count) noexcept
609{
610#if defined(DSPARK_SIMD_AVX)
611 int i = 0;
612 for (; i + 7 < count; i += 8)
613 {
614 __m256 vDst = _mm256_loadu_ps(dst + i);
615 __m256 vSrc = _mm256_loadu_ps(src + i);
616 _mm256_storeu_ps(dst + i, _mm256_add_ps(vDst, vSrc));
617 }
618 for (; i < count; ++i) dst[i] += src[i];
619
620#elif defined(DSPARK_SIMD_SSE2)
621 int i = 0;
622 for (; i + 3 < count; i += 4)
623 {
624 __m128 vDst = _mm_loadu_ps(dst + i);
625 __m128 vSrc = _mm_loadu_ps(src + i);
626 _mm_storeu_ps(dst + i, _mm_add_ps(vDst, vSrc));
627 }
628 for (; i < count; ++i) dst[i] += src[i];
629
630#elif defined(DSPARK_SIMD_NEON)
631 int i = 0;
632 for (; i + 3 < count; i += 4)
633 {
634 float32x4_t vDst = vld1q_f32(dst + i);
635 float32x4_t vSrc = vld1q_f32(src + i);
636 vst1q_f32(dst + i, vaddq_f32(vDst, vSrc));
637 }
638 for (; i < count; ++i) dst[i] += src[i];
639
640#else
641 for (int i = 0; i < count; ++i) dst[i] += src[i];
642#endif
643}
644
646inline void add(double* DSPARK_RESTRICT dst, const double* DSPARK_RESTRICT src, int count) noexcept
647{
648#if defined(DSPARK_SIMD_AVX)
649 int i = 0;
650 for (; i + 3 < count; i += 4)
651 _mm256_storeu_pd(dst + i, _mm256_add_pd(_mm256_loadu_pd(dst + i), _mm256_loadu_pd(src + i)));
652 for (; i < count; ++i) dst[i] += src[i];
653
654#elif defined(DSPARK_SIMD_SSE2)
655 int i = 0;
656 for (; i + 1 < count; i += 2)
657 {
658 __m128d vDst = _mm_loadu_pd(dst + i);
659 __m128d vSrc = _mm_loadu_pd(src + i);
660 _mm_storeu_pd(dst + i, _mm_add_pd(vDst, vSrc));
661 }
662 for (; i < count; ++i) dst[i] += src[i];
663
664#elif defined(DSPARK_SIMD_NEON)
665 int i = 0;
666 for (; i + 1 < count; i += 2)
667 vst1q_f64(dst + i, vaddq_f64(vld1q_f64(dst + i), vld1q_f64(src + i)));
668 for (; i < count; ++i) dst[i] += src[i];
669#else
670 for (int i = 0; i < count; ++i) dst[i] += src[i];
671#endif
672}
673
674// ============================================================================
675// multiply -- dst[i] = a[i] * b[i]
676// ============================================================================
677
679inline void multiply(float* DSPARK_RESTRICT dst, const float* DSPARK_RESTRICT a, const float* DSPARK_RESTRICT b, int count) noexcept
680{
681#if defined(DSPARK_SIMD_AVX)
682 int i = 0;
683 for (; i + 7 < count; i += 8)
684 _mm256_storeu_ps(dst + i, _mm256_mul_ps(_mm256_loadu_ps(a + i), _mm256_loadu_ps(b + i)));
685 for (; i < count; ++i) dst[i] = a[i] * b[i];
686#elif defined(DSPARK_SIMD_SSE2)
687 int i = 0;
688 for (; i + 3 < count; i += 4)
689 _mm_storeu_ps(dst + i, _mm_mul_ps(_mm_loadu_ps(a + i), _mm_loadu_ps(b + i)));
690 for (; i < count; ++i) dst[i] = a[i] * b[i];
691#elif defined(DSPARK_SIMD_NEON)
692 int i = 0;
693 for (; i + 3 < count; i += 4)
694 vst1q_f32(dst + i, vmulq_f32(vld1q_f32(a + i), vld1q_f32(b + i)));
695 for (; i < count; ++i) dst[i] = a[i] * b[i];
696#else
697 for (int i = 0; i < count; ++i) dst[i] = a[i] * b[i];
698#endif
699}
700
702inline void multiply(double* DSPARK_RESTRICT dst, const double* DSPARK_RESTRICT a, const double* DSPARK_RESTRICT b, int count) noexcept
703{
704#if defined(DSPARK_SIMD_AVX)
705 int i = 0;
706 for (; i + 3 < count; i += 4)
707 _mm256_storeu_pd(dst + i, _mm256_mul_pd(_mm256_loadu_pd(a + i), _mm256_loadu_pd(b + i)));
708 for (; i < count; ++i) dst[i] = a[i] * b[i];
709#elif defined(DSPARK_SIMD_SSE2)
710 int i = 0;
711 for (; i + 1 < count; i += 2)
712 _mm_storeu_pd(dst + i, _mm_mul_pd(_mm_loadu_pd(a + i), _mm_loadu_pd(b + i)));
713 for (; i < count; ++i) dst[i] = a[i] * b[i];
714
715#elif defined(DSPARK_SIMD_NEON)
716 int i = 0;
717 for (; i + 1 < count; i += 2)
718 vst1q_f64(dst + i, vmulq_f64(vld1q_f64(a + i), vld1q_f64(b + i)));
719 for (; i < count; ++i) dst[i] = a[i] * b[i];
720#else
721 for (int i = 0; i < count; ++i) dst[i] = a[i] * b[i];
722#endif
723}
724
725// ============================================================================
726// copyWithGain -- dst[i] = src[i] * gain
727// ============================================================================
728
730inline void copyWithGain(float* DSPARK_RESTRICT dst, const float* DSPARK_RESTRICT src, float gain, int count) noexcept
731{
732#if defined(DSPARK_SIMD_AVX)
733 const __m256 vGain = _mm256_set1_ps(gain);
734 int i = 0;
735 for (; i + 7 < count; i += 8)
736 _mm256_storeu_ps(dst + i, _mm256_mul_ps(_mm256_loadu_ps(src + i), vGain));
737 for (; i < count; ++i) dst[i] = src[i] * gain;
738#elif defined(DSPARK_SIMD_SSE2)
739 const __m128 vGain = _mm_set1_ps(gain);
740 int i = 0;
741 for (; i + 3 < count; i += 4)
742 _mm_storeu_ps(dst + i, _mm_mul_ps(_mm_loadu_ps(src + i), vGain));
743 for (; i < count; ++i) dst[i] = src[i] * gain;
744#elif defined(DSPARK_SIMD_NEON)
745 const float32x4_t vGain = vdupq_n_f32(gain);
746 int i = 0;
747 for (; i + 3 < count; i += 4)
748 vst1q_f32(dst + i, vmulq_f32(vld1q_f32(src + i), vGain));
749 for (; i < count; ++i) dst[i] = src[i] * gain;
750#else
751 for (int i = 0; i < count; ++i) dst[i] = src[i] * gain;
752#endif
753}
754
756inline void copyWithGain(double* DSPARK_RESTRICT dst, const double* DSPARK_RESTRICT src, double gain, int count) noexcept
757{
758#if defined(DSPARK_SIMD_AVX)
759 const __m256d vGain = _mm256_set1_pd(gain);
760 int i = 0;
761 for (; i + 3 < count; i += 4)
762 _mm256_storeu_pd(dst + i, _mm256_mul_pd(_mm256_loadu_pd(src + i), vGain));
763 for (; i < count; ++i) dst[i] = src[i] * gain;
764#elif defined(DSPARK_SIMD_SSE2)
765 const __m128d vGain = _mm_set1_pd(gain);
766 int i = 0;
767 for (; i + 1 < count; i += 2)
768 _mm_storeu_pd(dst + i, _mm_mul_pd(_mm_loadu_pd(src + i), vGain));
769 for (; i < count; ++i) dst[i] = src[i] * gain;
770
771#elif defined(DSPARK_SIMD_NEON)
772 const float64x2_t vGain = vdupq_n_f64(gain);
773 int i = 0;
774 for (; i + 1 < count; i += 2)
775 vst1q_f64(dst + i, vmulq_f64(vld1q_f64(src + i), vGain));
776 for (; i < count; ++i) dst[i] = src[i] * gain;
777#else
778 for (int i = 0; i < count; ++i) dst[i] = src[i] * gain;
779#endif
780}
781
782// ============================================================================
783// applyGainRamp -- data[i] *= gainStart + i * step (linear fade)
784// ============================================================================
785
794inline void applyGainRamp(float* DSPARK_RESTRICT data, float gainStart, float gainEnd, int count) noexcept
795{
796 if (count <= 0) return;
797 const float step = (gainEnd - gainStart) / static_cast<float>(count);
798#if defined(DSPARK_SIMD_AVX)
799 const __m256 vStep = _mm256_set1_ps(step * 8.0f);
800 __m256 vGain = _mm256_setr_ps(gainStart, gainStart + step,
801 gainStart + 2 * step, gainStart + 3 * step,
802 gainStart + 4 * step, gainStart + 5 * step,
803 gainStart + 6 * step, gainStart + 7 * step);
804 int i = 0;
805 for (; i + 7 < count; i += 8)
806 {
807 _mm256_storeu_ps(data + i, _mm256_mul_ps(_mm256_loadu_ps(data + i), vGain));
808 vGain = _mm256_add_ps(vGain, vStep);
809 }
810 for (; i < count; ++i) data[i] *= gainStart + step * static_cast<float>(i);
811
812#elif defined(DSPARK_SIMD_SSE2)
813 const __m128 vStep = _mm_set1_ps(step * 4.0f);
814 __m128 vGain = _mm_setr_ps(gainStart, gainStart + step, gainStart + 2 * step, gainStart + 3 * step);
815 int i = 0;
816 for (; i + 3 < count; i += 4)
817 {
818 _mm_storeu_ps(data + i, _mm_mul_ps(_mm_loadu_ps(data + i), vGain));
819 vGain = _mm_add_ps(vGain, vStep);
820 }
821 for (; i < count; ++i) data[i] *= gainStart + step * static_cast<float>(i);
822
823#elif defined(DSPARK_SIMD_NEON)
824 const float32x4_t vStep = vdupq_n_f32(step * 4.0f);
825 const float init[4] = { gainStart, gainStart + step, gainStart + 2 * step, gainStart + 3 * step };
826 float32x4_t vGain = vld1q_f32(init);
827 int i = 0;
828 for (; i + 3 < count; i += 4)
829 {
830 vst1q_f32(data + i, vmulq_f32(vld1q_f32(data + i), vGain));
831 vGain = vaddq_f32(vGain, vStep);
832 }
833 for (; i < count; ++i) data[i] *= gainStart + step * static_cast<float>(i);
834
835#else
836 for (int i = 0; i < count; ++i) data[i] *= gainStart + step * static_cast<float>(i);
837#endif
838}
839
841inline void applyGainRamp(double* DSPARK_RESTRICT data, double gainStart, double gainEnd, int count) noexcept
842{
843 if (count <= 0) return;
844 const double step = (gainEnd - gainStart) / static_cast<double>(count);
845#if defined(DSPARK_SIMD_AVX)
846 const __m256d vStep = _mm256_set1_pd(step * 4.0);
847 __m256d vGain = _mm256_setr_pd(gainStart, gainStart + step,
848 gainStart + 2 * step, gainStart + 3 * step);
849 int i = 0;
850 for (; i + 3 < count; i += 4)
851 {
852 _mm256_storeu_pd(data + i, _mm256_mul_pd(_mm256_loadu_pd(data + i), vGain));
853 vGain = _mm256_add_pd(vGain, vStep);
854 }
855 for (; i < count; ++i) data[i] *= gainStart + step * static_cast<double>(i);
856
857#elif defined(DSPARK_SIMD_SSE2)
858 const __m128d vStep = _mm_set1_pd(step * 2.0);
859 __m128d vGain = _mm_setr_pd(gainStart, gainStart + step);
860 int i = 0;
861 for (; i + 1 < count; i += 2)
862 {
863 _mm_storeu_pd(data + i, _mm_mul_pd(_mm_loadu_pd(data + i), vGain));
864 vGain = _mm_add_pd(vGain, vStep);
865 }
866 for (; i < count; ++i) data[i] *= gainStart + step * static_cast<double>(i);
867
868#elif defined(DSPARK_SIMD_NEON)
869 const float64x2_t vStep = vdupq_n_f64(step * 2.0);
870 const double init[2] = { gainStart, gainStart + step };
871 float64x2_t vGain = vld1q_f64(init);
872 int i = 0;
873 for (; i + 1 < count; i += 2)
874 {
875 vst1q_f64(data + i, vmulq_f64(vld1q_f64(data + i), vGain));
876 vGain = vaddq_f64(vGain, vStep);
877 }
878 for (; i < count; ++i) data[i] *= gainStart + step * static_cast<double>(i);
879
880#else
881 for (int i = 0; i < count; ++i) data[i] *= gainStart + step * static_cast<double>(i);
882#endif
883}
884
885// ============================================================================
886// addWithGainRamp -- dst[i] += src[i] * (gainStart + i * step)
887// ============================================================================
888
896inline void addWithGainRamp(float* DSPARK_RESTRICT dst, const float* DSPARK_RESTRICT src,
897 float gainStart, float gainEnd, int count) noexcept
898{
899 if (count <= 0) return;
900 const float step = (gainEnd - gainStart) / static_cast<float>(count);
901#if defined(DSPARK_SIMD_AVX)
902 const __m256 vStep = _mm256_set1_ps(step * 8.0f);
903 __m256 vGain = _mm256_setr_ps(gainStart, gainStart + step,
904 gainStart + 2 * step, gainStart + 3 * step,
905 gainStart + 4 * step, gainStart + 5 * step,
906 gainStart + 6 * step, gainStart + 7 * step);
907 int i = 0;
908 for (; i + 7 < count; i += 8)
909 {
910 const __m256 vDst = _mm256_loadu_ps(dst + i);
911 const __m256 vSrc = _mm256_loadu_ps(src + i);
912 #if defined(DSPARK_SIMD_FMA)
913 _mm256_storeu_ps(dst + i, _mm256_fmadd_ps(vSrc, vGain, vDst));
914 #else
915 _mm256_storeu_ps(dst + i, _mm256_add_ps(vDst, _mm256_mul_ps(vSrc, vGain)));
916 #endif
917 vGain = _mm256_add_ps(vGain, vStep);
918 }
919 for (; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<float>(i));
920
921#elif defined(DSPARK_SIMD_SSE2)
922 const __m128 vStep = _mm_set1_ps(step * 4.0f);
923 __m128 vGain = _mm_setr_ps(gainStart, gainStart + step, gainStart + 2 * step, gainStart + 3 * step);
924 int i = 0;
925 for (; i + 3 < count; i += 4)
926 {
927 const __m128 vDst = _mm_loadu_ps(dst + i);
928 const __m128 vSrc = _mm_loadu_ps(src + i);
929 #if defined(DSPARK_SIMD_FMA)
930 _mm_storeu_ps(dst + i, _mm_fmadd_ps(vSrc, vGain, vDst));
931 #else
932 _mm_storeu_ps(dst + i, _mm_add_ps(vDst, _mm_mul_ps(vSrc, vGain)));
933 #endif
934 vGain = _mm_add_ps(vGain, vStep);
935 }
936 for (; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<float>(i));
937
938#elif defined(DSPARK_SIMD_NEON)
939 const float32x4_t vStep = vdupq_n_f32(step * 4.0f);
940 const float init[4] = { gainStart, gainStart + step, gainStart + 2 * step, gainStart + 3 * step };
941 float32x4_t vGain = vld1q_f32(init);
942 int i = 0;
943 for (; i + 3 < count; i += 4)
944 {
945 vst1q_f32(dst + i, vfmaq_f32(vld1q_f32(dst + i), vld1q_f32(src + i), vGain));
946 vGain = vaddq_f32(vGain, vStep);
947 }
948 for (; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<float>(i));
949
950#else
951 for (int i = 0; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<float>(i));
952#endif
953}
954
956inline void addWithGainRamp(double* DSPARK_RESTRICT dst, const double* DSPARK_RESTRICT src,
957 double gainStart, double gainEnd, int count) noexcept
958{
959 if (count <= 0) return;
960 const double step = (gainEnd - gainStart) / static_cast<double>(count);
961#if defined(DSPARK_SIMD_AVX)
962 const __m256d vStep = _mm256_set1_pd(step * 4.0);
963 __m256d vGain = _mm256_setr_pd(gainStart, gainStart + step,
964 gainStart + 2 * step, gainStart + 3 * step);
965 int i = 0;
966 for (; i + 3 < count; i += 4)
967 {
968 const __m256d vDst = _mm256_loadu_pd(dst + i);
969 const __m256d vSrc = _mm256_loadu_pd(src + i);
970 #if defined(DSPARK_SIMD_FMA)
971 _mm256_storeu_pd(dst + i, _mm256_fmadd_pd(vSrc, vGain, vDst));
972 #else
973 _mm256_storeu_pd(dst + i, _mm256_add_pd(vDst, _mm256_mul_pd(vSrc, vGain)));
974 #endif
975 vGain = _mm256_add_pd(vGain, vStep);
976 }
977 for (; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<double>(i));
978
979#elif defined(DSPARK_SIMD_SSE2)
980 const __m128d vStep = _mm_set1_pd(step * 2.0);
981 __m128d vGain = _mm_setr_pd(gainStart, gainStart + step);
982 int i = 0;
983 for (; i + 1 < count; i += 2)
984 {
985 const __m128d vDst = _mm_loadu_pd(dst + i);
986 const __m128d vSrc = _mm_loadu_pd(src + i);
987 #if defined(DSPARK_SIMD_FMA)
988 _mm_storeu_pd(dst + i, _mm_fmadd_pd(vSrc, vGain, vDst));
989 #else
990 _mm_storeu_pd(dst + i, _mm_add_pd(vDst, _mm_mul_pd(vSrc, vGain)));
991 #endif
992 vGain = _mm_add_pd(vGain, vStep);
993 }
994 for (; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<double>(i));
995
996#elif defined(DSPARK_SIMD_NEON)
997 const float64x2_t vStep = vdupq_n_f64(step * 2.0);
998 const double init[2] = { gainStart, gainStart + step };
999 float64x2_t vGain = vld1q_f64(init);
1000 int i = 0;
1001 for (; i + 1 < count; i += 2)
1002 {
1003 vst1q_f64(dst + i, vfmaq_f64(vld1q_f64(dst + i), vld1q_f64(src + i), vGain));
1004 vGain = vaddq_f64(vGain, vStep);
1005 }
1006 for (; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<double>(i));
1007
1008#else
1009 for (int i = 0; i < count; ++i) dst[i] += src[i] * (gainStart + step * static_cast<double>(i));
1010#endif
1011}
1012
1013// ============================================================================
1014// sumOfSquares -- sum(data[i]^2) (RMS / energy metering)
1015// ============================================================================
1016
1018inline float sumOfSquares(const float* DSPARK_RESTRICT data, int count) noexcept
1019{
1020 return dotProduct(data, data, count);
1021}
1022
1024inline double sumOfSquares(const double* DSPARK_RESTRICT data, int count) noexcept
1025{
1026 return dotProduct(data, data, count);
1027}
1028
1029// ============================================================================
1030// complexMulAccum -- accum[k] += a[k] * b[k] over interleaved complex bins
1031// ============================================================================
1032
1039inline void complexMulAccum(float* DSPARK_RESTRICT accum, const float* DSPARK_RESTRICT a,
1040 const float* DSPARK_RESTRICT b, int bins) noexcept
1041{
1042#if defined(DSPARK_SIMD_AVX)
1043 // 4 complex bins per iteration. vpermilps builds the [re,re] / [im,im]
1044 // duplicates and the [im,re] swap within each 128-bit lane; the
1045 // alternating subtract/add of the two partial products comes straight
1046 // from (fm)addsub, so no sign mask is needed on this path.
1047 int k = 0;
1048 for (; k + 3 < bins; k += 4)
1049 {
1050 const __m256 va = _mm256_loadu_ps(a + 2 * k);
1051 const __m256 vb = _mm256_loadu_ps(b + 2 * k);
1052 const __m256 vacc = _mm256_loadu_ps(accum + 2 * k);
1053
1054 const __m256 aRe = _mm256_permute_ps(va, _MM_SHUFFLE(2, 2, 0, 0));
1055 const __m256 aIm = _mm256_permute_ps(va, _MM_SHUFFLE(3, 3, 1, 1));
1056 const __m256 bSwap = _mm256_permute_ps(vb, _MM_SHUFFLE(2, 3, 0, 1));
1057
1058 #if defined(DSPARK_SIMD_FMA)
1059 const __m256 prod = _mm256_fmaddsub_ps(aRe, vb, _mm256_mul_ps(aIm, bSwap));
1060 #else
1061 const __m256 prod = _mm256_addsub_ps(_mm256_mul_ps(aRe, vb), _mm256_mul_ps(aIm, bSwap));
1062 #endif
1063
1064 _mm256_storeu_ps(accum + 2 * k, _mm256_add_ps(vacc, prod));
1065 }
1066 for (; k < bins; ++k)
1067 {
1068 const float re1 = a[2 * k], im1 = a[2 * k + 1];
1069 const float re2 = b[2 * k], im2 = b[2 * k + 1];
1070 accum[2 * k] += re1 * re2 - im1 * im2;
1071 accum[2 * k + 1] += re1 * im2 + im1 * re2;
1072 }
1073#elif defined(DSPARK_SIMD_SSE2)
1074 // Sign mask: negate real positions for (re1*re2 - im1*im2)
1075 const __m128 negMask = _mm_castsi128_ps(_mm_setr_epi32(
1076 static_cast<int>(0x80000000u), 0,
1077 static_cast<int>(0x80000000u), 0));
1078
1079 int k = 0;
1080 for (; k + 1 < bins; k += 2)
1081 {
1082 const __m128 va = _mm_loadu_ps(a + 2 * k);
1083 const __m128 vb = _mm_loadu_ps(b + 2 * k);
1084 __m128 vacc = _mm_loadu_ps(accum + 2 * k);
1085
1086 const __m128 aRe = _mm_shuffle_ps(va, va, _MM_SHUFFLE(2, 2, 0, 0));
1087 const __m128 aIm = _mm_shuffle_ps(va, va, _MM_SHUFFLE(3, 3, 1, 1));
1088 const __m128 bSwap = _mm_shuffle_ps(vb, vb, _MM_SHUFFLE(2, 3, 0, 1));
1089
1090 const __m128 p1 = _mm_mul_ps(aRe, vb);
1091 const __m128 p2 = _mm_xor_ps(_mm_mul_ps(aIm, bSwap), negMask);
1092
1093 vacc = _mm_add_ps(vacc, _mm_add_ps(p1, p2));
1094 _mm_storeu_ps(accum + 2 * k, vacc);
1095 }
1096 for (; k < bins; ++k)
1097 {
1098 const float re1 = a[2 * k], im1 = a[2 * k + 1];
1099 const float re2 = b[2 * k], im2 = b[2 * k + 1];
1100 accum[2 * k] += re1 * re2 - im1 * im2;
1101 accum[2 * k + 1] += re1 * im2 + im1 * re2;
1102 }
1103#elif defined(DSPARK_SIMD_NEON)
1104 alignas(16) static constexpr uint32_t kNegRe[4] = { 0x80000000u, 0u, 0x80000000u, 0u };
1105 const uint32x4_t negMask = vld1q_u32(kNegRe);
1106
1107 int k = 0;
1108 for (; k + 1 < bins; k += 2)
1109 {
1110 const float32x4_t va = vld1q_f32(a + 2 * k);
1111 const float32x4_t vb = vld1q_f32(b + 2 * k);
1112 float32x4_t vacc = vld1q_f32(accum + 2 * k);
1113
1114 const float32x4_t aRe = vtrn1q_f32(va, va);
1115 const float32x4_t aIm = vtrn2q_f32(va, va);
1116 const float32x4_t bSwap = vrev64q_f32(vb);
1117
1118 const float32x4_t p1 = vmulq_f32(aRe, vb);
1119 const float32x4_t p2 = vreinterpretq_f32_u32(
1120 veorq_u32(vreinterpretq_u32_f32(vmulq_f32(aIm, bSwap)), negMask));
1121
1122 vacc = vaddq_f32(vacc, vaddq_f32(p1, p2));
1123 vst1q_f32(accum + 2 * k, vacc);
1124 }
1125 for (; k < bins; ++k)
1126 {
1127 const float re1 = a[2 * k], im1 = a[2 * k + 1];
1128 const float re2 = b[2 * k], im2 = b[2 * k + 1];
1129 accum[2 * k] += re1 * re2 - im1 * im2;
1130 accum[2 * k + 1] += re1 * im2 + im1 * re2;
1131 }
1132#else
1133 for (int k = 0; k < bins; ++k)
1134 {
1135 const float re1 = a[2 * k], im1 = a[2 * k + 1];
1136 const float re2 = b[2 * k], im2 = b[2 * k + 1];
1137 accum[2 * k] += re1 * re2 - im1 * im2;
1138 accum[2 * k + 1] += re1 * im2 + im1 * re2;
1139 }
1140#endif
1141}
1142
1144inline void complexMulAccum(double* DSPARK_RESTRICT accum, const double* DSPARK_RESTRICT a,
1145 const double* DSPARK_RESTRICT b, int bins) noexcept
1146{
1147 for (int k = 0; k < bins; ++k)
1148 {
1149 const double re1 = a[2 * k], im1 = a[2 * k + 1];
1150 const double re2 = b[2 * k], im2 = b[2 * k + 1];
1151 accum[2 * k] += re1 * re2 - im1 * im2;
1152 accum[2 * k + 1] += re1 * im2 + im1 * re2;
1153 }
1154}
1155
1156// ============================================================================
1157// Template dispatchers
1158// ============================================================================
1159
1160template <typename T>
1161void addWithGainT(T* DSPARK_RESTRICT dst, const T* DSPARK_RESTRICT src, T gain, int count) noexcept
1162{
1163 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1164 addWithGain(dst, src, gain, count);
1165}
1166
1167template <typename T>
1168void applyGainT(T* DSPARK_RESTRICT data, T gain, int count) noexcept
1169{
1170 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1171 applyGain(data, gain, count);
1172}
1173
1174template <typename T>
1175T peakLevelT(const T* DSPARK_RESTRICT data, int count) noexcept
1176{
1177 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1178 return peakLevel(data, count);
1179}
1180
1181template <typename T>
1182T dotProductT(const T* DSPARK_RESTRICT a, const T* DSPARK_RESTRICT b, int count) noexcept
1183{
1184 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1185 return dotProduct(a, b, count);
1186}
1187
1188template <typename T>
1189void addT(T* DSPARK_RESTRICT dst, const T* DSPARK_RESTRICT src, int count) noexcept
1190{
1191 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1192 add(dst, src, count);
1193}
1194
1195template <typename T>
1196void multiplyT(T* DSPARK_RESTRICT dst, const T* DSPARK_RESTRICT a, const T* DSPARK_RESTRICT b, int count) noexcept
1197{
1198 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1199 multiply(dst, a, b, count);
1200}
1201
1202template <typename T>
1203void copyWithGainT(T* DSPARK_RESTRICT dst, const T* DSPARK_RESTRICT src, T gain, int count) noexcept
1204{
1205 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1206 copyWithGain(dst, src, gain, count);
1207}
1208
1209template <typename T>
1210T sumOfSquaresT(const T* DSPARK_RESTRICT data, int count) noexcept
1211{
1212 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1213 return sumOfSquares(data, count);
1214}
1215
1216template <typename T>
1217void applyGainRampT(T* DSPARK_RESTRICT data, T gainStart, T gainEnd, int count) noexcept
1218{
1219 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1220 applyGainRamp(data, gainStart, gainEnd, count);
1221}
1222
1223template <typename T>
1224void addWithGainRampT(T* DSPARK_RESTRICT dst, const T* DSPARK_RESTRICT src, T gainStart, T gainEnd, int count) noexcept
1225{
1226 static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "SimdOps: only float and double are supported");
1227 addWithGainRamp(dst, src, gainStart, gainEnd, count);
1228}
1229
1230} // namespace simd
1231} // namespace dspark
DSPARK_RESTRICT
#define DSPARK_RESTRICT
Definition SimdOps.h:60
dspark::simd::peakLevelT
T peakLevelT(const T *DSPARK_RESTRICT data, int count) noexcept
Definition SimdOps.h:1175
dspark::simd::dotProductT
T dotProductT(const T *DSPARK_RESTRICT a, const T *DSPARK_RESTRICT b, int count) noexcept
Definition SimdOps.h:1182
dspark::simd::multiply
void multiply(float *DSPARK_RESTRICT dst, const float *DSPARK_RESTRICT a, const float *DSPARK_RESTRICT b, int count) noexcept
Element-wise product of two buffers into a destination.
Definition SimdOps.h:679
dspark::simd::complexMulAccum
void complexMulAccum(float *DSPARK_RESTRICT accum, const float *DSPARK_RESTRICT a, const float *DSPARK_RESTRICT b, int bins) noexcept
Complex multiply-accumulate over interleaved [re, im, ...] spectra.
Definition SimdOps.h:1039
dspark::simd::applyGainRamp
void applyGainRamp(float *DSPARK_RESTRICT data, float gainStart, float gainEnd, int count) noexcept
In-place linear gain ramp (the canonical click-free fade primitive).
Definition SimdOps.h:794
dspark::simd::addWithGain
void addWithGain(float *DSPARK_RESTRICT dst, const float *DSPARK_RESTRICT src, float gain, int count) noexcept
Adds source samples scaled by a gain factor into a destination buffer.
Definition SimdOps.h:77
dspark::simd::applyGainT
void applyGainT(T *DSPARK_RESTRICT data, T gain, int count) noexcept
Definition SimdOps.h:1168
dspark::simd::copyWithGainT
void copyWithGainT(T *DSPARK_RESTRICT dst, const T *DSPARK_RESTRICT src, T gain, int count) noexcept
Definition SimdOps.h:1203
dspark::simd::addT
void addT(T *DSPARK_RESTRICT dst, const T *DSPARK_RESTRICT src, int count) noexcept
Definition SimdOps.h:1189
dspark::simd::addWithGainRamp
void addWithGainRamp(float *DSPARK_RESTRICT dst, const float *DSPARK_RESTRICT src, float gainStart, float gainEnd, int count) noexcept
Adds source samples scaled by a linear gain ramp into a destination.
Definition SimdOps.h:896
dspark::simd::addWithGainRampT
void addWithGainRampT(T *DSPARK_RESTRICT dst, const T *DSPARK_RESTRICT src, T gainStart, T gainEnd, int count) noexcept
Definition SimdOps.h:1224
dspark::simd::sumOfSquaresT
T sumOfSquaresT(const T *DSPARK_RESTRICT data, int count) noexcept
Definition SimdOps.h:1210
dspark::simd::dotProduct
float dotProduct(const float *DSPARK_RESTRICT a, const float *DSPARK_RESTRICT b, int count) noexcept
Computes the dot product of two arrays.
Definition SimdOps.h:419
dspark::simd::applyGainRampT
void applyGainRampT(T *DSPARK_RESTRICT data, T gainStart, T gainEnd, int count) noexcept
Definition SimdOps.h:1217
dspark::simd::copyWithGain
void copyWithGain(float *DSPARK_RESTRICT dst, const float *DSPARK_RESTRICT src, float gain, int count) noexcept
Copies a buffer applying a gain factor (out-of-place).
Definition SimdOps.h:730
dspark::simd::add
void add(float *DSPARK_RESTRICT dst, const float *DSPARK_RESTRICT src, int count) noexcept
Adds source samples into a destination buffer (no scaling).
Definition SimdOps.h:608
dspark::simd::multiplyT
void multiplyT(T *DSPARK_RESTRICT dst, const T *DSPARK_RESTRICT a, const T *DSPARK_RESTRICT b, int count) noexcept
Definition SimdOps.h:1196
dspark::simd::addWithGainT
void addWithGainT(T *DSPARK_RESTRICT dst, const T *DSPARK_RESTRICT src, T gain, int count) noexcept
Definition SimdOps.h:1161
dspark::simd::sumOfSquares
float sumOfSquares(const float *DSPARK_RESTRICT data, int count) noexcept
Returns the sum of squared samples (energy).
Definition SimdOps.h:1018
dspark::simd::peakLevel
float peakLevel(const float *DSPARK_RESTRICT data, int count) noexcept
Returns the peak absolute sample value in a buffer.
Definition SimdOps.h:264
dspark::simd::applyGain
void applyGain(float *DSPARK_RESTRICT data, float gain, int count) noexcept
Multiplies all samples in a buffer by a gain factor.
Definition SimdOps.h:186
dspark
Main namespace for the DSPark framework.
Definition AnalogConstants.h:29