73# include <string_view>
74# include <type_traits>
76# if (defined(__GNUC__) or defined(__clang__)) \
77 and (defined(__x86_64__) or defined(__i386__) \
78 or defined(_M_X64) or defined(_M_IX86)) \
83# include <immintrin.h>
84# define ALEPH_FFT_HAS_X86_AVX2_DISPATCH 1
85# define ALEPH_FFT_AVX2_TARGET __attribute__((target("avx2")))
87# define ALEPH_FFT_HAS_X86_AVX2_DISPATCH 0
90# if (defined(__GNUC__) or defined(__clang__)) \
91 and (defined(__aarch64__) or defined(__arm__) or defined(_M_ARM64) \
93# if defined(__ARM_NEON) or defined(__ARM_NEON__)
95# define ALEPH_FFT_HAS_ARM_NEON_DISPATCH 1
97# define ALEPH_FFT_HAS_ARM_NEON_DISPATCH 0
100# define ALEPH_FFT_HAS_ARM_NEON_DISPATCH 0
103# if ALEPH_FFT_HAS_ARM_NEON_DISPATCH and defined(__aarch64__)
104# define ALEPH_FFT_HAS_ARM_NEON_DOUBLE 1
106# define ALEPH_FFT_HAS_ARM_NEON_DOUBLE 0
109# if ALEPH_FFT_HAS_ARM_NEON_DISPATCH and defined(__linux__)
110# include <sys/auxv.h>
111# include <asm/hwcap.h>
156 template <std::
floating_po
int Real =
double>
355 template <
typename T>
359 template <
typename T>
363 template <
typename Container>
371 template <
typename Container>
376 requires std::convertible_to<
decltype(*std::begin(c)),
Real>;
380 template <
typename Container>
385 requires std::convertible_to<std::remove_cvref_t<
decltype(*std::begin(c))>,
389 template <
typename Container>
403 return std::polar(
Real(1),
angle *
static_cast<Real>(index));
410 return std::max(
Real(1),
411 std::log2(
static_cast<Real>(std::max(n,
size_t(2)))));
422 * std::numeric_limits<Real>::epsilon() *
scale;
430 const size_t n)
noexcept
432 const Real scale = std::abs(lhs.real()) + std::abs(lhs.imag())
433 + std::abs(rhs.real()) + std::abs(rhs.imag())
436 * std::numeric_limits<Real>::epsilon() *
scale;
445 const size_t n =
input.size();
447 << ctx <<
": input must be non-empty";
451 << ctx <<
": coefficient 0 has non-negligible imaginary part "
457 const size_t half = n / 2;
461 << ctx <<
": coefficient " <<
half
462 <<
" has non-negligible imaginary part " <<
input[
half].imag();
470 << ctx <<
": spectrum is not Hermitian at k=" <<
k;
489 for (
size_t i = 0; i < n; ++i)
491 const Real phase =
Real(2) * std::numbers::pi_v<Real>
493 const Real value = a0 - a1 * std::cos(phase)
494 + a2 * std::cos(
Real(2) * phase);
502 template <
typename T>
508 for (
size_t i =
input.size(); i > 0; --i)
514 template <
typename T>
519 <<
"FFT::prefix_copy: length " << length
520 <<
" exceeds input size " <<
input.size();
524 for (
size_t i = 0; i < length; ++i)
530 template <
typename T>
537 <<
"FFT::slice_copy: offset " <<
offset <<
" exceeds input size " <<
input.size();
539 <<
"FFT::slice_copy: length " << length <<
" exceeds slice capacity " << (
input.size() -
offset);
543 for (
size_t i = 0; i < length; ++i)
558 > std::numeric_limits<size_t>::max() / 3 ?
559 std::numeric_limits<size_t>::max() :
574 <<
"FFT::reflect_pad_signal: signal size must be >= 2 when pad_len > 0";
576 <<
"FFT::reflect_pad_signal: pad length " <<
pad_len
577 <<
" must be smaller than signal size " << signal.
size();
582 for (
size_t i = 0; i <
pad_len; ++i)
585 for (
size_t i = 0; i < signal.
size(); ++i)
588 for (
size_t i = 0; i <
pad_len; ++i)
603 for (
size_t i = 0; i <
left_pad; ++i)
605 for (
size_t i = 0; i < signal.
size(); ++i)
619 const size_t hop_size,
623 << ctx <<
": analysis window must be non-empty";
625 << ctx <<
": synthesis window must be non-empty";
628 <<
" does not match synthesis window size "
633 for (
size_t i = 0; i < hop_size; ++i)
634 profile(i) =
Real(0);
649 *
Real(128) * std::numeric_limits<Real>::epsilon();
650 for (
size_t i = 0; i < profile.size(); ++i)
651 if (std::abs(profile[i]) <= tol)
662 if (profile.is_empty())
667 *
Real(256) * std::numeric_limits<Real>::epsilon();
668 for (
size_t i = 0; i < profile.size(); ++i)
669 if (std::abs(profile[i] -
mean) > tol)
678 const size_t hop_size,
679 const bool validate_nola,
680 const bool validate_cola,
683 if (
not validate_nola
and not validate_cola)
693 << ctx <<
": window pair does not satisfy NOLA for hop size "
698 << ctx <<
": window pair does not satisfy COLA for hop size "
705 const size_t frame_size,
706 const size_t hop_size,
736 for (
size_t i = 0; i <
input.size(); ++i)
746 for (
size_t i = 0; i <
input.size(); ++i)
756 for (
size_t i = 0; i <
input.size(); ++i)
765 if (
input.is_empty())
770 for (
size_t i = 0; i <
input.size(); ++i)
783 << ctx <<
": sample rate " <<
sample_rate <<
" must be strictly positive";
786 for (
size_t k = 0;
k < frequency.
size(); ++
k)
788 /
static_cast<Real>(fft_size);
795 const size_t fft_size)
noexcept
799 if (fft_size % 2 == 0
and bin == fft_size / 2)
807 const size_t frame_size,
810 const size_t hop_size =
811 options.hop_size == 0 ? std::max(
static_cast<size_t>(1), frame_size / 2) :
options.hop_size;
819 const size_t frame_size,
822 const size_t fft_size =
825 << ctx <<
": FFT size " << fft_size <<
" is smaller than frame size " << frame_size;
845 << ctx <<
": framing produced no analysis frames";
847 for (
size_t i = 0; i < frames.
size(); ++i)
858 template <
typename Density>
862 const size_t fft_size)
noexcept
872 if (
input.is_empty())
876 *
Real(64) * std::numeric_limits<Real>::epsilon();
877 size_t n =
input.size();
878 while (n > 1
and std::abs(
input[n - 1]) <= tol)
888 for (
size_t i = 0; i <
input.size(); ++i)
919 return (std::abs(reference) +
Real(1))
920 * multiplier * std::numeric_limits<Real>::epsilon();
929 * std::numeric_limits<Real>::epsilon();
941 << ctx <<
": numerator must be non-empty";
943 << ctx <<
": denominator must be non-empty";
945 const Real a0 = denominator[0];
947 *
Real(64) * std::numeric_limits<Real>::epsilon();
949 << ctx <<
": leading denominator coefficient must be non-zero";
956 for (
size_t i = 0; i <= order; ++i)
982 <<
" does not match a dense " << n <<
"x" << n <<
" system";
984 << ctx <<
": rhs size " << rhs.
size() <<
" does not match system size " << n;
987 for (
size_t i = 0; i <
matrix.size(); ++i)
989 for (
size_t i = 0; i < rhs.
size(); ++i)
993 *
Real(128) * std::numeric_limits<Real>::epsilon();
1015 << ctx <<
": singular or ill-conditioned linear system";
1019 for (
size_t k =
col;
k < n; ++
k)
1035 for (
size_t k =
col + 1;
k < n; ++
k)
1037 rhs(
row) -= factor * rhs[
col];
1042 for (
size_t i = 0; i < n; ++i)
1043 solution(i) =
Real(0);
1047 const size_t i =
row - 1;
1049 for (
size_t col = i + 1;
col < n; ++
col)
1054 << ctx <<
": singular or ill-conditioned linear system";
1067 << ctx <<
": numerator size " << numerator.
size()
1068 <<
" does not match denominator size " << denominator.
size();
1070 const size_t order = denominator.
size() - 1;
1075 for (
size_t i = 0; i <
system.size(); ++i)
1083 for (
size_t i = 0; i < order; ++i)
1086 for (
size_t i = 0; i < order; ++i)
1088 coeff(i, 0) += denominator[i + 1];
1094 for (
size_t i = 0; i < order; ++i)
1095 rhs(i) = numerator[i + 1] - denominator[i + 1] * numerator[0];
1109 << ctx <<
": numerator size " << numerator.
size()
1110 <<
" does not match denominator size " << denominator.
size();
1112 const size_t order = denominator.
size() - 1;
1116 <<
" does not match filter order " << order;
1124 for (
size_t i = 0; i < signal.
size(); ++i)
1125 output(i) = numerator[0] * signal[i];
1130 for (
size_t i = 0; i < order; ++i)
1133 for (
size_t n = 0; n < signal.
size(); ++n)
1135 const Real x = signal[n];
1136 const Real y = numerator[0] * x + state[0];
1139 for (
size_t i = 0; i + 1 < order; ++i)
1140 state(i) = state[i + 1] + numerator[i + 1] * x
1141 - denominator[i + 1] *
y;
1142 state(order - 1) = numerator[order] * x - denominator[order] *
y;
1167 for (
size_t i = 0; i < signal.
size(); ++i)
1168 output(i) = signal[i] * gain;
1206 for (
size_t i = 0; i <
sections.size(); ++i)
1223 for (
size_t j = 0; j <
output.size(); ++j)
1264 template <
typename T>
1269 for (
size_t i = 0; i < src.
size(); ++i)
1273 template <
typename T>
1278 <<
"FFT::drop_prefix: count " <<
count <<
" exceeds input size " <<
input.size();
1291 for (
size_t i =
count; i <
input.size(); ++i)
1299 return lhs == 0
or rhs == 0 ?
1301 lhs > std::numeric_limits<size_t>::max() / rhs ?
1302 std::numeric_limits<size_t>::max() :
1312 return pool !=
nullptr
1322 for (
size_t i = 0; i <
batch.size(); ++i)
1337 tile = std::max<size_t>(4, tile);
1338 tile = std::min(tile, std::max<size_t>(
static_cast<size_t>(1),
batch_size));
1339 if (pool !=
nullptr and pool->num_threads() > 1)
1340 tile = std::max(tile, pool->num_threads() * 2);
1341 return std::min(tile, std::max<size_t>(
static_cast<size_t>(1),
batch_size));
1346 const size_t frame_size)
1349 for (
size_t i = 0; i < frame_size; ++i)
1352 const size_t length = std::min(frame_size,
input.size());
1353 for (
size_t i = 0; i < length; ++i)
1354 frame(i) =
input[i];
1369 << ctx <<
": FFT size " << fft_size <<
" is smaller than window size " << window.
size();
1383 const size_t fft_size)
1391 const size_t fft_size,
1394 const size_t chunk_size)
1396 if (frames.is_empty())
1420 for (
size_t i = 0; i <
count; ++i)
1430 for (
size_t i = 0; i <
count; ++i)
1442 const size_t chunk_size = 0)
1447 const size_t fft_size =
1457 if (frames.is_empty())
1484 << ctx <<
": transfer response is singular at omega=" << omega;
1498 << ctx <<
": number of frequency samples must be positive";
1527 << ctx <<
": coefficient series must be non-empty";
1535 for (
size_t idx = n - 1; idx-- > 0;)
1548 [[
nodiscard]]
static PolynomialEvaluation
1558 for (
size_t i = 1; i < n; ++i)
1576 for (
size_t i = 1; i + 1 < n; ++i)
1625 std::numeric_limits<Real>::min());
1626 for (
size_t i = 1; i < monic.size(); ++i)
1627 bound = std::max(bound,
Real(1) + std::abs(monic[i]) /
leading);
1634 if (monic.size() <= 1)
1637 const Real constant = std::abs(monic[monic.size() - 1]);
1642 for (
size_t i = 0; i + 1 < monic.size(); ++i)
1654 << ctx <<
": polynomial degree must be at least one";
1660 << ctx <<
": leading polynomial coefficient must be non-zero";
1664 for (
size_t i = 0; i < n; ++i)
1675 if (
not std::isfinite(root_scale)
1676 or root_scale <= std::numeric_limits<Real>::min())
1677 root_scale =
Real(1);
1683 for (
size_t i = 1; i < n; ++i)
1689 output.root_scale = root_scale;
1714 for (
size_t i = 0; i < degree; ++i)
1717 Real(2) * std::numbers::pi_v<Real>
1718 * (
static_cast<Real>(i) +
Real(0.5))
1719 /
static_cast<Real>(degree);
1723 +
Real(0.2) *
static_cast<Real>(i + 1)
1724 /
static_cast<Real>(degree + 1));
1725 roots(i) = std::polar(radius,
angle);
1737 const size_t degree = roots.
size();
1741 for (
size_t i = 0; i < degree; ++i)
1746 "FFT::root_solver");
1748 for (
size_t j = 0; j < degree; ++j)
1751 const Complex delta = roots[i] - roots[j];
1752 if (std::abs(delta) > tol)
1757 if (std::abs(
denom) <= tol)
1762 roots(i) += std::polar(
nudge,
1764 * std::numbers::pi_v<Real>
1766 /
static_cast<Real>(degree + 1));
1775 if (std::abs(step) > tol * (
Real(1) + std::abs(roots[i])))
1792 const size_t degree = roots.
size();
1796 for (
size_t i = 0; i < degree; ++i)
1801 "FFT::root_solver");
1803 for (
size_t j = 0; j < degree; ++j)
1806 const Complex delta = roots[i] - roots[j];
1807 if (std::abs(delta) <= tol)
1810 roots(i) += std::polar(
nudge,
1812 * std::numbers::pi_v<Real>
1814 /
static_cast<Real>(degree + 1));
1821 if (std::abs(
denom) <= tol)
1829 if (std::abs(step) > tol * (
Real(1) + std::abs(roots[i])))
1844 const size_t iterations)
noexcept
1846 for (
size_t i = 0; i < roots.size(); ++i)
1852 "FFT::root_solver");
1858 if (std::abs(step) <= tol * (
Real(1) + std::abs(roots[i])))
1865 const Real tol)
noexcept
1868 for (
size_t i = 0; i < roots.size(); ++i)
1873 if (std::abs(roots[i].
imag()) <= tol * (
Real(8) + std::abs(roots[i])))
1880 const Complex target = std::conj(roots[i]);
1881 size_t best = roots.size();
1883 for (
size_t j = i + 1; j < roots.size(); ++j)
1887 const Real error = std::abs(roots[j] - target);
1895 if (
best == roots.size()
1900 const Real imag = std::abs(center.imag());
1901 roots(i) =
Complex(center.real(),
1903 roots(
best) = std::conj(roots[i]);
1912 const Real tol)
noexcept
1914 for (
size_t i = 0; i < roots.size(); ++i)
1919 "FFT::root_solver");
1921 const Real radius = std::abs(roots[i]);
1922 for (
size_t k = 0;
k < coefficients.size(); ++
k)
1923 bound = bound * radius + std::abs(coefficients[
k]);
1925 if (std::abs(eval.
value) > tol * (bound +
Real(1)))
1937 const size_t n = coefficients.size();
1938 const size_t degree = n - 1;
1941 for (
size_t i = 0; i < n; ++i)
1944 for (
size_t k = 0;
k < degree; ++
k)
1946 const size_t m = degree -
k;
1949 *
static_cast<Real>(
k)
1950 /
static_cast<Real>(degree + 1),
1951 Real(2) * std::numbers::pi_v<Real>
1953 /
static_cast<Real>(degree));
1961 for (
size_t i = 1; i <=
m; ++i)
1969 if (std::abs(p) <= tol)
1992 step = (
Real(1) + std::abs(z))
1993 * std::polar(
Real(1),
1994 Real(2) * std::numbers::pi_v<Real>
1999 if (std::abs(step) <= tol * (
Real(1) + std::abs(z)))
2013 for (
size_t i = 1; i <
m; ++i)
2038 << ctx <<
": leading polynomial coefficient must be non-zero";
2040 const size_t degree = n - 1;
2044 roots(0) =
Complex(-coefficients[1] / coefficients[0],
Real(0));
2050 const Real a = coefficients[0];
2051 const Real b = coefficients[1];
2052 const Real c = coefficients[2];
2059 if (std::abs(q) <= tol)
2097 for (
size_t i = 0; i < roots.
size(); ++i)
2098 roots(i) *=
problem.root_scale;
2107 << ctx <<
": root solver residual check failed";
2118 if (
input.is_empty())
2125 if (first ==
input.size())
2130 for (
size_t i = first; i <
input.size(); ++i)
2168 << ctx <<
": transfer response is singular at omega=" << omega;
2190 *
Real(256) * std::numeric_limits<Real>::epsilon();
2217 << ctx <<
": at least one biquad section is required";
2221 for (
size_t i = 0; i < response.
omega.
size(); ++i)
2225 for (
size_t j = 0; j <
sections.size(); ++j)
2265 *
Real(128) * std::numeric_limits<Real>::epsilon();
2278 for (
size_t i = 0; i < roots.size(); ++i)
2279 radius = std::max(radius,
static_cast<Real>(std::abs(roots[i])));
2302 for (
size_t i = 0; i <
zeros.
size(); ++i)
2307 for (
size_t j = 0; j <
poles.
size(); ++j)
2330 pair.has_zero =
true;
2331 pair.has_pole =
true;
2337 for (
size_t i = 0; i <
zeros.
size(); ++i)
2342 pair.has_zero =
true;
2346 for (
size_t i = 0; i <
poles.
size(); ++i)
2351 pair.has_pole =
true;
2364 Real best = std::numeric_limits<Real>::infinity();
2365 for (
size_t i = 0; i < pairs.size(); ++i)
2366 if (pairs[i].has_zero
and pairs[i].has_pole)
2367 best = std::min(
best, pairs[i].distance);
2405 template <
typename T>
2413 for (
size_t i = 0; i <
output.size(); ++i)
2416 for (
size_t i = 0; i < lhs.
size(); ++i)
2417 for (
size_t j = 0; j < rhs.
size(); ++j)
2418 output(i + j) += lhs[i] * rhs[j];
2424 template <
typename T>
2431 <<
"FFT::add_scaled_polynomial: destination size "
2432 <<
dst.size() <<
" is smaller than source size " << src.
size();
2434 for (
size_t i = 0; i < src.
size(); ++i)
2443 for (
size_t i = 0; i <
exponent; ++i)
2453 coeffs(0) =
Real(1);
2454 for (
size_t k = 1;
k <= power; ++
k)
2455 coeffs(
k) = coeffs[
k - 1]
2456 *
static_cast<Real>(power -
k + 1)
2457 /
static_cast<Real>(
k)
2472 << ctx <<
": analog polynomial must be non-empty";
2474 << ctx <<
": sample rate " <<
sample_rate <<
" must be positive";
2478 << ctx <<
": analog degree " << degree <<
" exceeds requested bilinear order " << order;
2481 for (
size_t i = 0; i <
output.size(); ++i)
2485 for (
size_t i = 0; i <
analog.size(); ++i)
2487 const size_t power = degree - i;
2506 for (
size_t i = 0; i < roots.
size(); ++i)
2509 for (
size_t j = 0; j <
next.size(); ++j)
2512 for (
size_t j = 0; j < coeffs.
size(); ++j)
2514 next(j) += coeffs[j];
2515 next(j + 1) -= coeffs[j] * roots[i];
2518 coeffs = std::move(
next);
2543 for (
size_t i = 0; i <
exponent; ++i)
2560 << ctx <<
": polynomial must be non-empty";
2562 << ctx <<
": rational substitution polynomials must be non-empty";
2564 const size_t order = degree - 1;
2569 for (
size_t i = 1; i <= order; ++i)
2578 for (
size_t i = 0; i <= order; ++i)
2586 for (
size_t j = 0; j <
grown.size(); ++j)
2588 for (
size_t j = 0; j <
output.size(); ++j)
2594 for (
size_t j = 0; j <
aligned.size(); ++j)
2596 for (
size_t j = 0; j <
term.size(); ++j)
2598 for (
size_t j = 0; j <
output.size(); ++j)
2649 << ctx <<
": transformed analog denominator leading coefficient "
2650 << a0 <<
" is too small";
2671 for (
size_t i = 0; i < coeffs.
size(); ++i)
2673 static_cast<Real>(std::abs(coeffs[i])));
2676 *
Real(1024) * std::numeric_limits<Real>::epsilon();
2678 for (
size_t i = 0; i < coeffs.
size(); ++i)
2681 << ctx <<
": root set does not generate a real polynomial";
2696 for (
size_t i = 0; i < roots.
size(); ++i)
2698 return center /
static_cast<Real>(roots.
size());
2709 *
Real(2048) * std::numeric_limits<Real>::epsilon();
2713 for (
size_t i = 0; i < result.
size(); ++i)
2720 const Complex target = std::conj(result[i]);
2721 for (
size_t j = i + 1; j < result.
size(); ++j)
2725 const Real error = std::abs(result[j] - target);
2736 (result[i].real() + result[
best].real()) /
Real(2);
2738 (std::abs(result[i].
imag())
2763 *
Real(2048) * std::numeric_limits<Real>::epsilon();
2767 for (
size_t i = 0; i < roots.
size(); ++i)
2769 if (
used[i]
or std::abs(roots[i].
imag()) <= tol)
2774 const Complex target = std::conj(roots[i]);
2775 for (
size_t j = i + 1; j < roots.
size(); ++j)
2780 const Real error = std::abs(roots[j] - target);
2790 + std::abs(roots[i])))
2791 << ctx <<
": complex root " << roots[i] <<
" does not have a matching conjugate";
2802 for (
size_t i = 0; i < roots.
size(); ++i)
2808 << ctx <<
": root " << roots[i] <<
" was not paired into a real section";
2816 for (
size_t i = 0; i <
real_roots.size(); ++i)
2828 for (
size_t j = first + 1; j <
real_roots.size(); ++j)
2867 << ctx <<
": section denominator order " << (coeffs.
denominator.
size() - 1)
2868 <<
" exceeds second order";
2873 for (
size_t i = 0; i < 3; ++i)
2901 << ctx <<
": analog numerator must be non-empty";
2903 << ctx <<
": analog denominator must be non-empty";
2905 const size_t order =
2931 << ctx <<
": sample rate " <<
sample_rate <<
" must be positive";
2934 <<
" must be positive";
2940 * std::tan(std::numbers::pi_v<Real>
2954# if defined(__GLIBCXX__) or defined(_MSVC_STL_VERSION)
2959 return std::ellint_1(
k,
phi);
2965 return std::comp_ellint_1(
k);
2976 for (
size_t i = 0; i < 64
and std::abs(a - b) >
2977 std::numeric_limits<Real>::epsilon() * a; ++i)
2980 b = std::sqrt(a * b);
2983 return std::numbers::pi_v<Real> / (
Real(2) * a);
2993 std::pow(std::numeric_limits<Real>::epsilon() /
Real(3),
2995 for (
size_t i = 0; i < 64; ++i)
3001 if (std::max({std::abs(dx), std::abs(dy), std::abs(
dz)})
3011 const Real lam = std::sqrt(x *
y) + std::sqrt(x * z)
3017 return Real(1) / std::sqrt((x +
y + z) /
Real(3));
3028 std::numeric_limits<Real>::epsilon() *
Real(4))
3031 if (std::abs(
phi) <=
3032 std::numeric_limits<Real>::epsilon() *
Real(4))
3039 Real(1) -
k *
k * s * s),
3053 << ctx <<
": elliptic modulus " <<
modulus <<
" must lie in [0, 1)";
3066 << ctx <<
": elliptic modulus " <<
modulus <<
" must lie in [0, 1)";
3076 values.
sn = std::sin(u);
3077 values.
cn = std::cos(u);
3080 values.
sn = -values.
sn;
3091 values.
dn = std::sqrt(std::max(
Real(0),
3094 values.
sn = -values.
sn;
3099 << ctx <<
": Jacobi argument " << u
3100 <<
" exceeds the supported quarter period " << K;
3106 std::numeric_limits<Real>::epsilon(),
3107 upper - std::numeric_limits<Real>::epsilon());
3116 if (residual >
Real(0))
3122 const Real derivative =
3123 Real(1) / std::sqrt(std::max(
Real(0),
3135 values.
sn = std::sin(
phi);
3136 values.
cn = std::cos(
phi);
3137 values.
dn = std::sqrt(std::max(
Real(0),
3140 * values.
sn * values.
sn));
3142 values.
sn = -values.
sn;
3152 << ctx <<
": inverse Jacobi sc value " <<
value
3153 <<
" must be non-negative";
3155 << ctx <<
": elliptic modulus " <<
modulus
3156 <<
" must lie in [0, 1)";
3166 << ctx <<
": order must be positive";
3168 << ctx <<
": selectivity modulus " <<
k1
3169 <<
" must lie in (0, 1)";
3175 / (
static_cast<Real>(order)
3179 Real upper =
Real(1) - std::sqrt(std::numeric_limits<Real>::epsilon());
3212 << ctx <<
": elliptic cd denominator vanished";
3221 << ctx <<
": order must be positive";
3226 for (
size_t k = 0;
k < order; ++
k)
3229 std::numbers::pi_v<Real>
3232 +
static_cast<Real>(order))
3233 / (
Real(2) *
static_cast<Real>(order));
3246 << ctx <<
": order must be positive";
3248 << ctx <<
": passband ripple " <<
ripple_db
3249 <<
" must be positive in decibels";
3251 const Real epsilon =
3253 const Real mu = std::asinh(
Real(1) / epsilon)
3254 /
static_cast<Real>(order);
3258 for (
size_t k = 0;
k < order; ++
k)
3261 std::numbers::pi_v<Real>
3264 / (
Real(2) *
static_cast<Real>(order));
3267 std::cosh(
mu) * std::cos(
theta));
3281 const Real attenuation_db,
3285 << ctx <<
": order must be positive";
3287 << ctx <<
": stopband attenuation " << attenuation_db
3288 <<
" must be positive in decibels";
3290 const Real epsilon =
3291 Real(1) / std::sqrt(std::pow(
Real(10), attenuation_db /
Real(10))
3293 const Real mu = std::asinh(
Real(1) / epsilon)
3294 /
static_cast<Real>(order);
3299 for (
size_t k = 0;
k < order; ++
k)
3302 std::numbers::pi_v<Real>
3304 / (
Real(2) *
static_cast<Real>(order));
3306 std::cosh(
mu) * std::cos(
theta));
3310 for (
size_t k = 0;
k < order / 2; ++
k)
3313 std::numbers::pi_v<Real>
3315 / (
Real(2) *
static_cast<Real>(order));
3327 denominator[denominator.
size() - 1] / numerator[numerator.
size() - 1];
3334 const Real attenuation_db,
3338 << ctx <<
": order must be positive";
3340 << ctx <<
": passband ripple " <<
ripple_db
3341 <<
" must be positive in decibels";
3343 << ctx <<
": stopband attenuation " << attenuation_db
3344 <<
" must exceed passband ripple " <<
ripple_db;
3346 const Real epsilon =
3349 std::pow(
Real(10), attenuation_db /
Real(10)) -
Real(1);
3351 << ctx <<
": stopband attenuation produced a non-positive stop term";
3371 K *
r / (
static_cast<Real>(order) *
K1);
3382 /
static_cast<Real>(order);
3387 << ctx <<
": elliptic prototype zero scale vanished";
3395 if (pole.real() >
Real(0))
3397 if (pole.imag() <
Real(0))
3398 pole = std::conj(pole);
3401 prototype.poles.append(std::conj(pole));
3404 if ((order & 1) != 0)
3411 << ctx <<
": elliptic prototype real pole denominator vanished";
3422 Real(1) / std::sqrt(
Real(1) + epsilon * epsilon) :
3426 / numerator[numerator.
size() - 1];
3448 << ctx <<
": analog transfer response is singular at omega="
3461 << ctx <<
": order must be positive";
3464 for (
size_t k = 0;
k <= order; ++
k)
3466 const long double numerator =
3467 std::tgamma(
static_cast<long double>(2 * order -
k + 1));
3468 const long double denominator =
3469 std::pow(2.0L,
static_cast<long double>(order -
k))
3470 * std::tgamma(
static_cast<long double>(order -
k + 1))
3471 * std::tgamma(
static_cast<long double>(
k + 1));
3473 static_cast<Real>(numerator / denominator);
3497 << ctx <<
": could not bracket the -3 dB normalization point";
3504 const Real magnitude =
3509 if (magnitude > target)
3517 for (
size_t i = 0; i <
poles.
size(); ++i)
3554 << ctx <<
": prototype order " <<
prototype.poles.size()
3555 <<
" does not match requested order " << order;
3568 for (
size_t i = 0; i < order; ++i)
3579 const size_t degree =
pole_groups[i].roots.size();
3605 [[
nodiscard]]
static std::pair<Real, Real>
3613 <<
" must be positive";
3684 << ctx <<
": could not assign finite-zero group " << i
3685 <<
" to any SOS section";
3714 << ctx <<
": specialized band-stop path expects prototypes without zeros";
3729 {
Real(1),
Real(0), center * center},
3744 for (
size_t i = 0; i < order; ++i)
3797 for (
size_t k = 1;
k < 80; ++
k)
3800 / (
static_cast<Real>(
k) *
static_cast<Real>(
k));
3802 if (std::abs(
term) <= std::numeric_limits<Real>::epsilon() *
sum)
3812 const Real pix = std::numbers::pi_v<Real> * x;
3813 if (std::abs(
pix) <= std::numeric_limits<Real>::epsilon())
3815 return std::sin(
pix) /
pix;
3839 << ctx <<
": FIR gain vanished at normalization omega=" << omega;
3840 for (
size_t i = 0; i < coeffs.
size(); ++i)
3853 << ctx <<
": num_taps must be positive";
3855 << ctx <<
": window size " << window.
size()
3856 <<
" does not match num_taps " <<
num_taps;
3858 << ctx <<
": sample rate " <<
sample_rate <<
" must be positive";
3861 <<
" must be positive";
3869 for (
size_t n = 0; n <
num_taps; ++n)
3871 const Real m =
static_cast<Real>(n) - center;
3875 coeffs(n) =
ideal * window[n];
3917 return Real(0.5) * width * width;
3933 if (width <=
Real(0))
3953 << ctx <<
": num_taps must be positive";
3955 << ctx <<
": num_taps must be odd for a Type-I linear-phase design";
3957 << ctx <<
": sample rate " <<
sample_rate <<
" must be positive";
3959 << ctx <<
": bands must contain an even number of band edges";
3961 << ctx <<
": desired size " << desired.
size()
3962 <<
" does not match band edge count " <<
bands.size();
3966 << ctx <<
": weights size " << weights.
size()
3967 <<
" does not match band count " <<
band_count;
3973 << ctx <<
": first band edge " <<
bands[0] <<
" must start at 0";
3975 << ctx <<
": last band edge " <<
bands[
bands.size() - 1]
3976 <<
" must reach Nyquist " <<
nyquist;
3987 << ctx <<
": weight " << weights[i]
3988 <<
" at band " << i <<
" must be positive";
3994 for (
size_t i = 0; i <
bands.size(); ++i)
3998 << ctx <<
": band edge " << edge <<
" lies outside [0, Nyquist]";
4000 << ctx <<
": band edges must be non-decreasing";
4003 << ctx <<
": band " << (i / 2)
4004 <<
" has non-positive width";
4010 const size_t order = (
num_taps - 1) / 2;
4011 const size_t dimension = order + 1;
4014 for (
size_t i = 0; i <
system.size(); ++i)
4016 for (
size_t i = 0; i < rhs.
size(); ++i)
4027 for (
size_t row = 0;
row < dimension; ++
row)
4051 for (
size_t i = 0; i < coeffs.
size(); ++i)
4052 coeffs(i) =
Real(0);
4054 const size_t center = order;
4059 coeffs(center -
k) =
tap;
4060 coeffs(center +
k) =
tap;
4066 [[
nodiscard]]
static WeightedFrequencyGrid
4076 << ctx <<
": num_taps must be positive";
4078 << ctx <<
": num_taps must be odd for a Type-I linear-phase design";
4080 << ctx <<
": sample rate " <<
sample_rate <<
" must be positive";
4083 <<
" must be at least 8";
4085 << ctx <<
": bands must contain an even number of band edges";
4087 << ctx <<
": desired size " << desired.
size()
4088 <<
" does not match band edge count " <<
bands.size();
4092 << ctx <<
": weights size " << weights.
size()
4093 <<
" does not match band count " <<
band_count;
4098 << ctx <<
": first band edge " <<
bands[0] <<
" must start at 0";
4100 << ctx <<
": last band edge " <<
bands[
bands.size() - 1]
4101 <<
" must reach Nyquist " <<
nyquist;
4113 << ctx <<
": band edges must lie in [0, Nyquist]";
4115 << ctx <<
": band edges must be non-decreasing";
4117 << ctx <<
": band " <<
band <<
" has non-positive width";
4123 << ctx <<
": weight " << weight <<
" at band " <<
band
4124 <<
" must be positive";
4134 / std::numbers::pi_v<Real>));
4145 static_cast<Real>(i)
4150 grid.omega.append(omega);
4152 grid.weight.append(weight);
4159 << ctx <<
": weighted frequency grid is too small";
4165 const Real omega)
noexcept
4168 for (
size_t k = 0;
k < coefficients.size(); ++
k)
4169 value += coefficients[
k] * std::cos(
static_cast<Real>(
k) * omega);
4181 const Real position =
4184 static_cast<Real>(i)
4187 extrema.append(
static_cast<size_t>(std::llround(position)));
4238 for (
size_t i = 1; i <
candidates.size(); ++i)
4240 if (
candidates[i].
sign == compressed[compressed.size() - 1].sign)
4242 if (
candidates[i].magnitude > compressed[compressed.size() - 1].magnitude)
4243 compressed(compressed.size() - 1) =
candidates[i];
4255 for (
size_t i = 0; i < compressed.size(); ++i)
4256 extrema.append(compressed[i].index);
4263 for (
size_t start = 0; start +
required_count <= compressed.size(); ++start)
4265 Real floor = std::numeric_limits<Real>::infinity();
4269 floor = std::min(
floor, compressed[start + i].magnitude);
4270 sum += compressed[start + i].magnitude;
4301 << ctx <<
": expected " << dimension <<
" extrema but got "
4306 for (
size_t i = 0; i <
matrix.size(); ++i)
4308 for (
size_t i = 0; i < rhs.
size(); ++i)
4311 for (
size_t row = 0;
row < dimension; ++
row)
4316 std::cos(
static_cast<Real>(
col) * omega[index]);
4318 ((
row & 1) == 0 ?
Real(1) :
Real(-1)) / weight[index];
4319 rhs(
row) = desired[index];
4336 << ctx <<
": max_iterations must be positive";
4352 for (
size_t i = 0; i <
cosine.size(); ++i)
4364 for (
size_t i = 0; i <
cosine.size(); ++i)
4368 for (
size_t i = 0; i <
grid.omega.size(); ++i)
4389 for (
size_t i = 0; i < coeffs.
size(); ++i)
4390 coeffs(i) =
Real(0);
4393 coeffs(center) =
cosine[0];
4394 for (
size_t k = 1;
k <
cosine.size(); ++
k)
4397 coeffs(center -
k) =
tap;
4398 coeffs(center +
k) =
tap;
4409 const Real target)
noexcept
4412 if (std::abs(dy) <= std::numeric_limits<Real>::epsilon())
4414 const Real t = (target -
y0) / dy;
4415 return x0 + t * (x1 -
x0);
4423 const Real x)
noexcept
4426 if (std::abs(dx) <= std::numeric_limits<Real>::epsilon())
4428 const Real t = (x -
x0) / dx;
4429 return y0 + t * (
y1 -
y0);
4438 / (
Real(2) * std::numbers::pi_v<Real>));
4442 template <
typename Evaluator>
4456 if (std::abs(
left_y) <= std::numeric_limits<Real>::epsilon())
4458 if (std::abs(
right_y) <= std::numeric_limits<Real>::epsilon())
4470 if (std::abs(
mid_y) <= std::numeric_limits<Real>::epsilon()
4496 << ctx <<
": at least one biquad section is required";
4499 for (
size_t i = 0; i <
sections.size(); ++i)
4511 <<
"FFT::phase_margin: omega size " << response.
omega.
size()
4512 <<
" does not match response size " << response.
response.
size();
4524 for (
size_t i = 1; i < response.
omega.
size(); ++i)
4526 const Real left = magnitude[i - 1] -
Real(1);
4527 const Real right = magnitude[i] -
Real(1);
4567 <<
"FFT::gain_margin: omega size " << response.
omega.
size()
4568 <<
" does not match response size " << response.
response.
size();
4580 for (
size_t i = 1; i < response.
omega.
size(); ++i)
4582 const Real p0 = phase[i - 1];
4583 const Real p1 = phase[i];
4587 static_cast<long long>(std::ceil((-
upper - std::numbers::pi_v<Real>)
4589 * std::numbers::pi_v<Real>)));
4591 static_cast<long long>(std::floor((-
lower - std::numbers::pi_v<Real>)
4593 * std::numbers::pi_v<Real>)));
4599 * std::numbers::pi_v<Real>;
4610 const Real ratio = std::abs(
mag) <= std::numeric_limits<Real>::epsilon() ?
4611 std::numeric_limits<Real>::infinity() :
4626 std::numeric_limits<Real>::infinity() :
4636 template <
typename Evaluator>
4642 <<
"FFT::phase_margin: omega size " << response.
omega.
size()
4643 <<
" does not match response size " << response.
response.
size();
4655 for (
size_t i = 1; i < response.
omega.
size(); ++i)
4657 const Real left = magnitude[i - 1] -
Real(1);
4658 const Real right = magnitude[i] -
Real(1);
4668 response.
omega[i - 1],
4702 template <
typename Evaluator>
4708 <<
"FFT::gain_margin: omega size " << response.
omega.
size()
4709 <<
" does not match response size " << response.
response.
size();
4721 for (
size_t i = 1; i < response.
omega.
size(); ++i)
4723 const Real p0 = phase[i - 1];
4724 const Real p1 = phase[i];
4728 static_cast<long long>(std::ceil((-
upper - std::numbers::pi_v<Real>)
4730 * std::numbers::pi_v<Real>)));
4732 static_cast<long long>(std::floor((-
lower - std::numbers::pi_v<Real>)
4734 * std::numbers::pi_v<Real>)));
4740 * std::numbers::pi_v<Real>;
4744 const Real reference =
4753 response.
omega[i - 1],
4760 std::abs(
mag) <= std::numeric_limits<Real>::epsilon() ?
4761 std::numeric_limits<Real>::infinity() :
4776 std::numeric_limits<Real>::infinity() :
4795 for (
size_t i = 1; i < phase.
size(); ++i)
4797 const Real delta = phase[i] - phase[i - 1];
4798 if (delta > std::numbers::pi_v<Real>)
4799 offset -=
Real(2) * std::numbers::pi_v<Real>;
4800 else if (delta < -std::numbers::pi_v<Real>)
4801 offset +=
Real(2) * std::numbers::pi_v<Real>;
4812 <<
"FFT::group_delay: omega size " << response.
omega.
size()
4813 <<
" does not match response size " << response.
response.
size();
4826 for (
size_t i = 0; i < response.
omega.
size(); ++i)
4828 const size_t left = i == 0 ? 0 : i - 1;
4829 const size_t right = i + 1 == response.
omega.
size() ? i : i + 1;
4831 if (std::abs(
domega) <= std::numeric_limits<Real>::epsilon())
4834 delay(i) = -(phase[right] - phase[left]) /
domega;
4845 <<
"FFT::phase_delay: omega size " << response.
omega.
size()
4846 <<
" does not match response size " << response.
response.
size();
4855 *
Real(128) * std::numeric_limits<Real>::epsilon();
4857 for (
size_t i = 0; i < response.
omega.
size(); ++i)
4869 const size_t fft_size,
4874 << ctx <<
": analysis window must be non-empty";
4876 << ctx <<
": synthesis window must be non-empty";
4879 <<
" does not match synthesis window size "
4882 << ctx <<
": hop size must be positive";
4884 << ctx <<
": FFT size must be positive";
4886 << ctx <<
": frame FFT size " << fft_size
4904 const size_t chunk_size)
4918 <<
"FFT::istft: frame " << i <<
" has size "
4919 <<
spectrogram[i].size() <<
" but expected " << fft_size;
4931 <<
" exceeds overlap-add length " <<
raw_length
4932 <<
" after centered trimming";
4958 *
Real(256) * std::numeric_limits<Real>::epsilon();
4965 <<
"FFT::istft: overlap-add normalization vanished at sample " << i;
4982 <<
"FFT::zero_padded_copy: target size " << n
4983 <<
" is smaller than input size " <<
input.
size();
4987 for (
size_t i = 0; i <
input.size(); ++i)
4989 for (
size_t i =
input.size(); i < n; ++i)
4999 <<
"FFT::zero_padded_copy: target size " << n
5000 <<
" is smaller than input size " <<
input.
size();
5004 for (
size_t i = 0; i <
input.size(); ++i)
5006 for (
size_t i =
input.size(); i < n; ++i)
5016 <<
"FFT::trim_to_size: target size " << n
5017 <<
" exceeds input size " <<
input.size();
5019 while (
input.size() > n)
5020 static_cast<void>(
input.remove_last());
5023 template <
typename T,
typename Container>
5028 requires std::convertible_to<
decltype(*std::begin(c)),
T>;
5035 if constexpr (
requires {
input.size(); })
5038 for (
auto it = std::begin(
input); it != std::end(
input); ++it)
5039 output.append(
static_cast<T>(*it));
5044 template <
typename Container>
5052 template <
typename Container>
5065 << ctx <<
": tensor shape must be non-empty";
5068 for (
size_t i = 0; i < shape.
size(); ++i)
5071 << ctx <<
": shape dimension " << i <<
" must be positive";
5084 for (
size_t i = shape.
size(); i > 0; --i)
5086 strides(i - 1) = stride;
5098 << ctx <<
": shape rank " << shape.
size()
5099 <<
" does not match strides rank " << strides.
size();
5101 size_t max_offset = 0;
5102 for (
size_t i = 0; i < shape.
size(); ++i)
5105 << ctx <<
": shape dimension " << i <<
" must be positive";
5123 << ctx <<
": tensor shape must be non-empty";
5126 if (
layout.strides.is_empty())
5130 << ctx <<
": flat data size " << data.
size()
5131 <<
" does not match row-major tensor size " <<
logical_size;
5137 << ctx <<
": strides rank " <<
normalized.strides.size()
5138 <<
" does not match shape rank " <<
normalized.shape.size();
5140 const size_t max_offset =
5143 << ctx <<
": tensor layout touches offset " << max_offset
5144 <<
" but flat buffer size is " << data.
size();
5157 << ctx <<
": at least one axis is required";
5162 for (
size_t i = 0; i <
axes.size(); ++i)
5165 << ctx <<
": axis " <<
axes[i] <<
" is out of range for rank "
5168 << ctx <<
": axis " <<
axes[i] <<
" appears more than once";
5187 for (
size_t i = 0; i < shape.
size(); ++i)
5203 for (
size_t i = 0; i <
counters.size(); ++i)
5210 for (
size_t i = 0; i <
other_dims.size(); ++i)
5214 for (
size_t rev =
other_dims.size(); rev > 0; --rev)
5216 const size_t idx = rev - 1;
5254 for (
size_t i = 0; i < slice.
size(); ++i)
5267 const size_t chunk_size = 0)
5274 "FFT::transform_axis");
5282 if (offsets.
size() == 1)
5286 if (pool !=
nullptr and pool->num_threads() > 1)
5287 plan.ptransform(*pool, slice, invert, chunk_size);
5289 plan.transform(slice, invert);
5299 plan.transform(slice, invert);
5304 and pool->num_threads() > 1
5308 for (
size_t i = 0; i < offsets.
size(); ++i)
5321 const size_t chunk_size = 0)
5349 for (
size_t i = 1; i <
rows; ++i)
5351 << ctx <<
": row " << i <<
" has size " <<
input[i].size()
5352 <<
" but expected " <<
cols;
5365 for (
size_t row = 0;
row < shape[0]; ++
row)
5366 for (
size_t col = 0;
col < shape[1]; ++
col)
5379 << ctx <<
": matrix shape must be positive";
5381 << ctx <<
": flat size " <<
input.size()
5382 <<
" does not match matrix shape " <<
rows <<
"x" <<
cols;
5402 << ctx <<
": tensor must be non-empty";
5406 << ctx <<
": tensor middle dimension must be positive";
5409 << ctx <<
": tensor innermost dimension must be positive";
5411 for (
size_t i = 0; i <
dim0; ++i)
5414 << ctx <<
": slab " << i <<
" has size " <<
input[i].size()
5415 <<
" but expected " <<
dim1;
5416 for (
size_t j = 0; j <
dim1; ++j)
5418 << ctx <<
": slab " << i <<
", row " << j
5419 <<
" has size " <<
input[i][j].size()
5420 <<
" but expected " <<
dim2;
5434 for (
size_t i = 0; i < shape[0]; ++i)
5435 for (
size_t j = 0; j < shape[1]; ++j)
5436 for (
size_t k = 0;
k < shape[2]; ++
k)
5450 << ctx <<
": tensor shape must be positive";
5452 << ctx <<
": flat size " <<
input.size()
5453 <<
" does not match tensor shape";
5457 for (
size_t i = 0; i <
dim0; ++i)
5460 for (
size_t j = 0; j <
dim1; ++j)
5463 for (
size_t k = 0;
k <
dim2; ++
k, ++index)
5479 if (source == target
or input.is_empty())
5485 const size_t frames =
input[0].size();
5487 frames == 0 ? 0 :
input[0][0].size();
5492 << ctx <<
": channel " <<
ch <<
" has " <<
input[
ch].size()
5493 <<
" frames but expected " << frames;
5494 for (
size_t frame = 0; frame < frames; ++frame)
5496 << ctx <<
": channel " <<
ch <<
", frame " << frame
5497 <<
" has " <<
input[
ch][frame].size()
5498 <<
" bins but expected " <<
bins;
5502 for (
size_t frame = 0; frame < frames; ++frame)
5514 for (
size_t frame = 0; frame < frames; ++frame)
5517 << ctx <<
": frame " << frame <<
" has " <<
input[frame].size()
5518 <<
" channels but expected " <<
channels;
5521 << ctx <<
": frame " << frame <<
", channel " <<
ch
5522 <<
" has " <<
input[frame][
ch].size()
5523 <<
" bins but expected " <<
bins;
5530 for (
size_t frame = 0; frame < frames; ++frame)
5557 const size_t chunk_size = 0)
5561 <<
"FFT::transform: input must be non-empty";
5573 const size_t half = n / 2;
5576 for (
size_t i = 0; i <
half; ++i)
5584 /
static_cast<Real>(n);
5597 if (pool !=
nullptr and pool->num_threads() > 1
and half > 1)
5600 for (
size_t k = 0;
k <
half; ++
k)
5610 template <
typename InversePackedTransform>
5615 const size_t chunk_size,
5620 << ctx <<
": input must be non-empty";
5625 return {
input[0].real()};
5629 << ctx <<
": input size " << n <<
" must be a power of two";
5633 const size_t half = n / 2;
5636 /
static_cast<Real>(n);
5651 for (
size_t k = 0;
k <
half; ++
k)
5666 for (
size_t i = 0; i <
half; ++i)
5679 const size_t chunk_size = 0)
5685 <<
"FFT::multiply: product size exceeds size_t capacity";
5694 for (
size_t i = 0; i < n; ++i)
5703 const size_t mirror = (n -
k) % n;
5711 if (pool !=
nullptr and pool->num_threads() > 1
and n > 1)
5714 for (
size_t k = 0;
k < n; ++
k)
5726 const size_t chunk_size = 0)
5732 <<
"FFT::multiply: product size exceeds size_t capacity";
5748 if (pool !=
nullptr and pool->num_threads() > 1
and n > 1)
5751 for (
size_t i = 0; i < n; ++i)
5778 for (
size_t i = 0; i <
input.size(); ++i)
5783 * std::numeric_limits<Real>::epsilon() *
scale;
5786 << ctx <<
": coefficient " << i
5787 <<
" has non-negligible imaginary part " <<
input[i].imag();
5806 const size_t n = a.size();
5807 for (
size_t i = 1, j = 0; i < n; ++i)
5809 size_t bit = n >> 1;
5810 for (; j & bit; bit >>= 1)
5814 std::swap(a(i), a(j));
5830 <<
"FFT::next_power_of_two: input size must be positive";
5836 <<
"FFT::next_power_of_two: size overflow for input " << n;
5847 const size_t chunk_size = 0)
5849 const size_t n = a.
size();
5851 <<
"FFT::transform: input must be non-empty";
5857 <<
"FFT::transform: input size " << n <<
" must be a power of two";
5862 for (
size_t tmp = n;
tmp > 1;
tmp >>= 1)
5871 const size_t blocks = n >> 1;
5872 auto r2_block = [&a](
const size_t block)
5874 const size_t base = block << 1;
5876 const Complex v = a[base + 1];
5878 a(base + 1) = u - v;
5881 if (pool !=
nullptr and pool->num_threads() > 1
and blocks > 1)
5884 for (
size_t block = 0; block < blocks; ++block)
5896 const size_t quarter = len >> 2;
5897 const size_t blocks = n / len;
5899 * std::numbers::pi_v<Real>
5900 /
static_cast<Real>(len);
5904 (
const size_t block)
5906 const size_t base = block * (
quarter << 2);
5908 for (
size_t j = 0; j <
quarter; ++j)
5913 const size_t i0 = base + j;
5938 const size_t next_j = j + 1;
5947 if (pool !=
nullptr and pool->num_threads() > 1
and blocks > 1)
5950 for (
size_t block = 0; block < blocks; ++block)
5957 for (
size_t i = 0; i < n; ++i)
5974 else if (n % 5 == 0)
5976 else if (n % 3 == 0)
5978 else if (n % 2 == 0)
6002 return n > 0
and (n & (n - 1)) == 0;
6006 [[
nodiscard]]
static constexpr const char *
6022 [[
nodiscard]]
static constexpr const char *
6043# if ALEPH_FFT_HAS_X86_AVX2_DISPATCH
6044 return std::same_as<Real, double>;
6054# if ALEPH_FFT_HAS_ARM_NEON_DISPATCH
6055 if constexpr (std::same_as<Real, float>)
6057# if ALEPH_FFT_HAS_ARM_NEON_DOUBLE
6058 if constexpr (std::same_as<Real, double>)
6071# if ALEPH_FFT_HAS_X86_AVX2_DISPATCH
6072 if constexpr (std::same_as<Real, double>)
6074 static const bool available = []()
noexcept
6091# if ALEPH_FFT_HAS_ARM_NEON_DISPATCH
6092 if constexpr (std::same_as<Real, float>)
6097# if defined(__linux__)
6098 static const bool available = []()
noexcept
6103# if defined(__aarch64__) and defined(HWCAP_ASIMD)
6105# elif defined(__arm__) and defined(HWCAP_NEON)
6117# if ALEPH_FFT_HAS_ARM_NEON_DOUBLE
6118 if constexpr (std::same_as<Real, double>)
6123# if defined(__linux__) and defined(HWCAP_ASIMD)
6124 static const bool available = []()
noexcept
6169 const char *
disabled = std::getenv(
"ALEPH_FFT_DISABLE_AVX2");
6173 if (
const char *
mode = std::getenv(
"ALEPH_FFT_SIMD");
6177 if (
value ==
"scalar")
6179 if (
value ==
"avx2")
6181 if (
value ==
"neon")
6186 const char *enabled = std::getenv(
"ALEPH_FFT_ENABLE_AVX2");
6187 if (enabled !=
nullptr and enabled[0] !=
'\0' and enabled[0] !=
'0')
6328 for (
size_t i = 1, j = 0; i <
n_; ++i)
6330 size_t bit =
n_ >> 1;
6331 for (; j & bit; bit >>= 1)
6340 const size_t half =
static_cast<size_t>(1) <<
stage;
6341 const size_t len =
half << 1;
6343 /
static_cast<Real>(len);
6345 for (
size_t j = 0; j <
half; ++j)
6358 for (
size_t i = 0; i <
n_; ++i)
6360 Real(-2) * std::numbers::pi_v<Real>
6361 *
static_cast<Real>(i)
6362 /
static_cast<Real>(
n_));
6373 for (
size_t i = 0; i <
n_; ++i)
6391 for (
size_t i = 1; i <
n_; ++i)
6408 for (
size_t i = 0; i <
n_; ++i)
6413# if ALEPH_FFT_HAS_X86_AVX2_DISPATCH
6417 static_assert(
sizeof(
Complex) ==
sizeof(
double) * 2);
6423 const __m256d value)
requires std::same_as<Real, double>
6425 static_assert(
sizeof(
Complex) ==
sizeof(
double) * 2);
6438 const __m256d rhs)
requires std::same_as<Real, double>
6449 const bool invert)
requires std::same_as<Real, double>
6459 requires std::same_as<Real, double>
6471 requires std::same_as<Real, double>
6473 for (
size_t block = 0;
block < blocks; ++
block)
6483 const size_t block)
const
6484 requires std::same_as<Real, double>
6488 for (; j + 1 <
quarter; j += 2)
6499 const size_t i0 = base + j;
6529 const size_t i0 = base + j;
6555 const Real inv_n)
requires std::same_as<Real, double>
6559 for (; i + 1 < a.size(); i += 2)
6562 for (; i < a.size(); ++i)
6569 const size_t blocks,
6572 const bool invert)
const
6573 requires std::same_as<Real, double>
6575 for (
size_t block = 0;
block < blocks; ++
block)
6580# if ALEPH_FFT_HAS_ARM_NEON_DISPATCH
6587# if ALEPH_FFT_HAS_ARM_NEON_DOUBLE
6598 static_assert(
sizeof(
Complex) ==
sizeof(
float) * 2);
6600 vld2_f32(
reinterpret_cast<const float *
>(ptr));
6608 static_assert(
sizeof(
Complex) ==
sizeof(
float) * 2);
6645 const bool invert)
requires std::same_as<Real, float>
6652# if ALEPH_FFT_HAS_ARM_NEON_DOUBLE
6656 static_assert(
sizeof(
Complex) ==
sizeof(
double) * 2);
6658 vld2q_f64(
reinterpret_cast<const double *
>(ptr));
6666 static_assert(
sizeof(
Complex) ==
sizeof(
double) * 2);
6703 const bool invert)
requires std::same_as<Real, double>
6717 const size_t block)
const
6721 for (; j + 1 <
quarter; j += 2)
6732 const size_t i0 = base + j;
6762 const size_t i0 = base + j;
6790 if constexpr (std::same_as<Real, float>)
6793 for (; i + 1 < a.
size(); i += 2)
6801# if ALEPH_FFT_HAS_ARM_NEON_DOUBLE
6802 else if constexpr (std::same_as<Real, double>)
6805 for (; i + 1 < a.
size(); i += 2)
6814 for (; i<a.
size();++i)
6821 const size_t blocks,
6824 const bool invert)
const
6826 for (
size_t block = 0;
block < blocks; ++
block)
6844 if constexpr (
not std::same_as<Real, double>)
6850 if (pool !=
nullptr and pool->num_threads() > 1)
6875 if constexpr (
not (std::same_as<Real, double>
or std::same_as<Real, float>))
6881 if (pool !=
nullptr and pool->num_threads() > 1)
6902 const size_t chunk_size,
6908# if ALEPH_FFT_HAS_X86_AVX2_DISPATCH
6909 if constexpr (std::same_as<Real, double>)
6924 const size_t quarter =
static_cast<size_t>(1) <<
s_low;
6925 const size_t blocks =
n_ >> (
s_high + 1);
6939# if ALEPH_FFT_HAS_ARM_NEON_DISPATCH
6946 if constexpr (std::same_as<Real, double>
or std::same_as<Real, float>)
6953 const size_t blocks =
n_ >> 1;
6954 for (
size_t block = 0; block < blocks; ++block)
6956 const size_t base = block << 1;
6958 const Complex v = a[base + 1];
6960 a(base + 1) = u - v;
6970 const size_t blocks =
n_ >> (
s_high + 1);
6989 const size_t blocks =
n_ >> 1;
6990 auto r2_block = [&a](
const size_t block)
6992 const size_t base = block << 1;
6994 const Complex v = a[base + 1];
6996 a(base + 1) = u - v;
7002 for (
size_t block = 0; block < blocks; ++block)
7015 const size_t quarter =
static_cast<size_t>(1) <<
s_low;
7016 const size_t blocks =
n_ >> (
s_high + 1);
7021 (
const size_t block)
7023 const size_t base = block * (
quarter << 2);
7024 for (
size_t j = 0; j <
quarter; ++j)
7032 const size_t i0 = base + j;
7059 for (
size_t block = 0; block < blocks; ++block)
7068 for (
size_t i = 0; i <
n_; ++i)
7077 const size_t chunk_size,
7081 <<
"FFT::Plan::transform: input size " << a.
size()
7082 <<
" does not match plan size " <<
n_;
7107 const bool invert)
const
7112 const size_t stride =
n_ / length;
7113 const size_t index = ((
exponent % length) * stride) %
n_;
7115 return invert ? std::conj(
root) :
root;
7121 const size_t current_size,
7122 const bool invert)
const
7124 if (current_size <= 1)
7128 const size_t sub_size = current_size / factor;
7132 for (
size_t q = 0; q < factor; ++q)
7136 for (
size_t t = 0; t <
sub_size; ++t)
7146 for (
size_t p = 0; p < factor; ++p)
7150 for (
size_t q = 0; q < factor; ++q)
7166 for (
size_t i = 0; i <
n_; ++i)
7176 const size_t chunk_size)
const
7180 auto initialize = [
this, &a, &work, invert](
const size_t i)
7189 work(i) = a[i] *
chirp;
7198 if (pool !=
nullptr)
7205 auto pointwise = [&work, &kernel](
const size_t i)
7207 work(i) *= kernel[i];
7216 if (pool !=
nullptr)
7221 auto finalize = [
this, &a, &work, invert](
const size_t i)
7224 a(i) = work[i] *
chirp;
7226 a(i) /=
static_cast<Real>(
n_);
7232 for (
size_t i = 0; i <
n_; ++i)
7259 <<
"FFT::Plan: size must be positive";
7305 const bool invert =
false)
const
7328 <<
"FFT::Plan::inverse_transform_real: input size " <<
input.
size()
7329 <<
" does not match plan size " <<
n_;
7336 input,
"FFT::Plan::inverse_transform_real",
nullptr, 0,
7351 const size_t chunk_size = 0)
const
7354 <<
"FFT::Plan::pinverse_transform_real: input size "
7355 <<
input.
size() <<
" does not match plan size " <<
n_;
7362 input,
"FFT::Plan::inverse_transform_real", &pool,
7389 const size_t chunk_size = 0)
const
7397 const bool invert =
false,
7398 const size_t chunk_size = 0)
const
7408 const size_t chunk_size = 0)
const
7421 for (
size_t i = 0; i <
batch.size(); ++i)
7424 <<
"FFT::Plan::transform_batch: batch item " << i
7425 <<
" has size " <<
batch[i].size()
7426 <<
" but plan size is " <<
n_;
7441 const size_t chunk_size = 0,
7446 for (
size_t i = 0; i <
batch.size(); ++i)
7448 <<
"FFT::Plan::ptransform_batch: batch item " << i
7449 <<
" has size " <<
batch[i].size()
7450 <<
" but plan size is " <<
n_;
7452 if (
batch.is_empty())
7457 for (
size_t i = 0; i <
batch.size(); ++i)
7476 const bool invert =
false,
7488 const bool invert =
false,
7489 const size_t chunk_size = 0,
7508 const size_t chunk_size = 0,
7521 for (
size_t i = 0; i <
input.size(); ++i)
7524 <<
"FFT::Plan::inverse_transform_real_batch: batch item " << i
7525 <<
" has size " <<
input[i].size()
7526 <<
" but plan size is " <<
n_;
7538 const size_t chunk_size = 0,
7541 for (
size_t i = 0; i <
input.size(); ++i)
7543 <<
"FFT::Plan::pinverse_transform_real_batch: batch item " << i
7544 <<
" has size " <<
input[i].size()
7545 <<
" but plan size is " <<
n_;
7547 if (
input.is_empty())
7553 for (
size_t i = 0; i <
input.size(); ++i)
7577 <<
"FFT::Plan::rfft: input size " <<
input.
size()
7578 <<
" does not match plan size " <<
n_;
7585 const size_t chunk_size = 0)
const
7588 <<
"FFT::Plan::prfft: input size " <<
input.
size()
7589 <<
" does not match plan size " <<
n_;
7606 "FFT::Plan::irfft"));
7612 const size_t chunk_size = 0)
const
7617 "FFT::Plan::pirfft"),
7633 const size_t chunk_size = 0)
const
7644 spectra,
n_,
"FFT::Plan::irfft_batch"));
7651 const size_t chunk_size = 0)
const
7657 "FFT::Plan::pirfft_batch"),
7667 const size_t chunk_size = 0)
7670 if (pool !=
nullptr)
7671 plan.ptransform(*pool,
lifted,
false, chunk_size);
7681 const size_t chunk_size = 0)
7685 for (
size_t i = 0; i <
input.size(); ++i)
7688 if (pool !=
nullptr)
7689 plan.ptransform_batch(*pool,
lifted,
false, chunk_size);
7699 <<
"FFT::compact_real_spectrum: spectrum must be non-empty";
7724 << ctx <<
": compact spectrum must be non-empty";
7729 << ctx <<
": compact spectrum size " <<
spectrum.
size()
7738 > std::numeric_limits<size_t>::max() / 2)
7739 << ctx <<
": inferred signal size overflows size_t";
7749 const size_t half = n / 2;
7752 for (
size_t i = 0; i < n; ++i)
7755 for (
size_t i = 0; i <=
half; ++i)
7757 for (
size_t i =
half + 1; i < n; ++i)
7758 full(i) = std::conj(full[n - i]);
7770 for (
size_t i = 0; i <
spectra.size(); ++i)
7781 for (
size_t i = 0; i <
input.size(); ++i)
7790 const size_t chunk_size = 0)
7792 const size_t n = a.
size();
7805 if (pool !=
nullptr)
7806 plan.ptransform(*pool, a, invert, chunk_size);
7808 plan.transform(a, invert);
7815 const size_t chunk_size)
7858 const size_t chunk_size = 0)
7896 const bool invert =
false,
const size_t chunk_size = 0)
7907 if (
batch.is_empty())
7919 const size_t chunk_size = 0)
7921 if (
batch.is_empty())
7925 plan.ptransform_batch(pool,
batch, invert, chunk_size);
7931 const bool invert =
false)
7942 const bool invert =
false,
7943 const size_t chunk_size = 0)
7961 const size_t chunk_size = 0)
7970 if (
input.is_empty())
7981 const size_t chunk_size = 0)
7983 if (
input.is_empty())
7987 return plan.prfft_batch(pool,
input, chunk_size);
8007 const size_t chunk_size = 0)
8013 return plan.pirfft_batch(pool,
spectra, chunk_size);
8043 const size_t chunk_size = 0)
8065 const size_t chunk_size = 0)
8075 const bool invert =
false)
8088 const bool invert =
false,
8089 const size_t chunk_size = 0)
8099 const bool invert =
false)
8107 "FFT::transformed2d");
8114 const bool invert =
false,
8115 const size_t chunk_size = 0)
8128 "FFT::ptransformed2d");
8142 const size_t chunk_size = 0)
8150 const bool invert =
false)
8154 "FFT::transformed2d_batch");
8160 "FFT::transformed2d_batch");
8167 const bool invert =
false,
8168 const size_t chunk_size = 0)
8172 "FFT::ptransformed2d_batch");
8183 "FFT::ptransformed2d_batch");
8189 const bool invert =
false)
8201 "FFT::transformed3d");
8208 const bool invert =
false,
8209 const size_t chunk_size = 0)
8223 "FFT::ptransformed3d");
8237 const size_t chunk_size = 0)
8253 "FFT::transpose_spectrogram_layout");
8283 const size_t chunk_size = 0)
8329 const size_t chunk_size = 0)
8356 const size_t chunk_size = 0)
8377 template <
typename Container>
8395 template <
typename Container>
8399 const size_t chunk_size = 0)
8417 template <
typename Container>
8436 template <
typename Container>
8440 const bool invert =
false,
const size_t chunk_size = 0)
8474 const size_t chunk_size = 0)
8494 const size_t chunk_size = 0)
8509 const size_t chunk_size = 0)
8516 template <
typename Container>
8525 template <
typename Container>
8529 const size_t chunk_size = 0)
8535 template <
typename Container>
8544 template <
typename Container>
8548 const size_t chunk_size = 0)
8554 template <
typename Container>
8563 template <
typename Container>
8567 const size_t chunk_size = 0)
8573 template <
typename Container>
8582 template <
typename Container>
8586 const size_t chunk_size = 0)
8601 const size_t chunk_size = 0)
8616 const size_t chunk_size = 0)
8622 template <
typename Container>
8631 template <
typename Container>
8635 const size_t chunk_size = 0)
8641 template <
typename Container>
8650 template <
typename Container>
8654 const size_t chunk_size = 0)
8678 const size_t chunk_size = 0)
8684 template <
typename Container>
8693 template <
typename Container>
8697 const size_t chunk_size = 0)
8718 "FFT::inverse_transform_real",
nullptr, 0,
8725 "FFT::inverse_transform_real",
8733 const size_t chunk_size = 0)
8737 "FFT::inverse_transform_real", &pool, chunk_size,
8744 "FFT::inverse_transform_real",
8768 const size_t chunk_size = 0)
8778 template <
typename Container>
8787 template <
typename Container>
8791 const size_t chunk_size = 0)
8797 template <
typename Container>
8806 template <
typename Container>
8812 const size_t chunk_size = 0)
8826 for (
size_t i = 0; i <
input.size(); ++i)
8832 template <
typename Container>
8846 for (
size_t i = 0; i <
input.size(); ++i)
8851 template <
typename Container>
8865 for (
size_t i = 0; i <
input.size(); ++i)
8870 template <
typename Container>
8904 <<
"FFT::kaiser_beta: attenuation " << attenuation_db
8905 <<
" must be non-negative";
8907 if (attenuation_db >
Real(50))
8908 return Real(0.1102) * (attenuation_db -
Real(8.7));
8909 if (attenuation_db >=
Real(21))
8910 return Real(0.5842) * std::pow(attenuation_db -
Real(21),
Real(0.4))
8911 +
Real(0.07886) * (attenuation_db -
Real(21));
8921 <<
"FFT::kaiser_window: beta " << beta <<
" must be non-negative";
8932 for (
size_t i = 0; i < n; ++i)
8935 const Real arg = beta * std::sqrt(std::max(
Real(0),
Real(1) - x * x));
8946 <<
"FFT::apply_window: signal size " << signal.
size()
8947 <<
" does not match window size " << window.
size();
8950 for (
size_t i = 0; i < signal.
size(); ++i)
8951 output(i) = signal[i] * window[i];
8960 <<
"FFT::apply_window: signal size " << signal.
size()
8961 <<
" does not match window size " << window.
size();
8964 for (
size_t i = 0; i < signal.
size(); ++i)
8965 output(i) = signal[i] * window[i];
8969 template <
typename SignalContainer,
typename WindowContainer>
8977 template <
typename SignalContainer,
typename WindowContainer>
9038 "FFT::firwin_lowpass");
9046 const Real attenuation_db)
9066 "FFT::firwin_highpass");
9067 const size_t center = (
num_taps - 1) / 2;
9068 for (
size_t i = 0; i < coeffs.
size(); ++i)
9069 coeffs(i) = -coeffs[i];
9070 coeffs(center) +=
Real(1);
9072 std::numbers::pi_v<Real>,
9073 "FFT::firwin_highpass");
9082 const Real attenuation_db)
9100 <<
" must be positive";
9113 "FFT::firwin_bandpass");
9119 "FFT::firwin_bandpass");
9121 for (
size_t i = 0; i <
num_taps; ++i)
9122 coeffs(i) =
low[i] - high[i];
9128 "FFT::firwin_bandpass");
9138 const Real attenuation_db)
9161 const size_t center = (
num_taps - 1) / 2;
9162 for (
size_t i = 0; i < coeffs.
size(); ++i)
9163 coeffs(i) = -coeffs[i];
9164 coeffs(center) +=
Real(1);
9175 const Real attenuation_db)
9200 template <
typename BandContainer,
typename DesiredContainer>
9255 template <
typename BandContainer,
typename DesiredContainer>
9298 const size_t up = 1,
9299 const size_t down = 1)
9302 <<
"FFT::upfirdn: up factor must be positive";
9304 <<
"FFT::upfirdn: down factor must be positive";
9306 <<
"FFT::upfirdn: FIR coefficients must be non-empty";
9311 <<
"FFT::upfirdn: upsampled signal length overflows size_t";
9319 for (
size_t n = 0; n < signal.
size(); ++n)
9321 const size_t base = n *
up;
9322 for (
size_t k = 0;
k < coeffs.
size(); ++
k)
9332 template <
typename SignalContainer,
typename CoeffContainer>
9337 const size_t up = 1,
9338 const size_t down = 1)
9354 <<
"FFT::resample_poly: up factor must be positive";
9356 <<
"FFT::resample_poly: down factor must be positive";
9358 <<
"FFT::resample_poly: FIR coefficients must be non-empty";
9362 const size_t gcd = std::gcd(
up,
down);
9369 const size_t expected_size =
9373 <<
"FFT::resample_poly: internal trim offset exceeds raw output";
9375 const size_t available = raw.
size() - start;
9376 const size_t trimmed_size = std::min(expected_size, available);
9378 if (
output.size() == expected_size)
9382 for (
size_t i = 0; i < expected_size; ++i)
9395 <<
"FFT::resample_poly: taps_per_phase must be positive";
9397 <<
"FFT::resample_poly: attenuation " <<
options.attenuation_db
9398 <<
" must be strictly positive";
9402 const size_t gcd = std::gcd(
up,
down);
9415 for (
size_t i = 0; i < coeffs.
size(); ++i)
9421 template <
typename SignalContainer,
typename CoeffContainer>
9435 template <
typename SignalContainer>
9446 template <
typename Container>
9454 template <
typename Container>
9462 template <
typename Container>
9470 template <
typename Container>
9478 template <
typename Container>
9486 template <
typename Container>
9508 template <
typename SignalContainer,
typename WindowContainer>
9517 template <
typename SignalContainer,
typename WindowContainer>
9539 <<
"FFT::window_coherent_gain: window must be non-empty";
9548 <<
"FFT::window_enbw: window must be non-empty";
9551 <<
"FFT::window_enbw: coherent gain is zero";
9556 template <
typename Container>
9564 template <
typename Container>
9572 template <
typename Container>
9583 const size_t frame_size,
9584 const size_t hop_size,
9585 const bool pad_end =
true)
9591 "FFT::frame_offsets");
9597 const size_t hop_size,
9598 const size_t signal_length = 0)
9601 <<
"FFT::overlap_add_frames: hop size must be positive";
9602 if (frames.is_empty())
9605 const size_t frame_size = frames[0].
size();
9607 <<
"FFT::overlap_add_frames: frames must be non-empty";
9608 for (
size_t i = 1; i < frames.size(); ++i)
9610 <<
"FFT::overlap_add_frames: frame " << i
9611 <<
" has size " << frames[i].size()
9612 <<
" but expected " << frame_size;
9614 const size_t total_length = hop_size * (frames.size() - 1) + frame_size;
9616 <<
"FFT::overlap_add_frames: requested signal length "
9617 << signal_length <<
" exceeds overlap-add extent " <<
total_length;
9623 for (
size_t frame = 0; frame < frames.size(); ++frame)
9625 const size_t offset = frame * hop_size;
9626 for (
size_t i = 0; i < frame_size; ++i)
9630 return signal_length == 0 ?
9636 [[
nodiscard]]
static PowerSpectralDensity
9642 const char *ctx =
"FFT::welch";
9647 << ctx <<
": window energy must be strictly positive";
9649 Plan
plan(fft_size);
9651 for (
size_t k = 0;
k < density.size(); ++
k)
9652 density(
k) =
Real(0);
9654 for (
size_t i = 0; i < frames.
size(); ++i)
9669 *
static_cast<Real>(frames.
size()));
9677 [[
nodiscard]]
static PowerSpectralDensity
9679 const size_t frame_size,
9687 [[
nodiscard]]
static CrossSpectralDensity
9694 const auto ctx =
"FFT::csd";
9696 << ctx <<
": signal sizes " << x.
size()
9697 <<
" and " <<
y.size() <<
" do not match";
9703 << ctx <<
": frame count mismatch after preparation";
9707 << ctx <<
": window energy must be strictly positive";
9709 Plan
plan(fft_size);
9711 for (
size_t k = 0;
k < density.size(); ++
k)
9714 for (
size_t i = 0; i <
x_frames.size(); ++i)
9718 for (
size_t k = 0;
k < density.size(); ++
k)
9737 [[
nodiscard]]
static CrossSpectralDensity
9740 const size_t frame_size,
9760 or pxy.density.size() !=
pyy.density.size())
9761 <<
"FFT::coherence: inconsistent spectrum sizes";
9763 CoherenceEstimate
output;
9766 const Real tol =
Real(1024) * std::numeric_limits<Real>::epsilon();
9767 for (
size_t k = 0;
k <
pxy.density.size(); ++
k)
9780 const size_t frame_size,
9787 template <
typename SignalContainer,
typename WindowContainer>
9789 [[
nodiscard]]
static PowerSpectralDensity
9801 template <
typename ContainerX,
typename ContainerY,
typename WindowContainer>
9805 [[
nodiscard]]
static CrossSpectralDensity
9819 template <
typename ContainerX,
typename ContainerY,
typename WindowContainer>
9841 const size_t hop_size)
9846 "FFT::window_overlap_profile");
9853 const size_t hop_size)
9864 const size_t hop_size)
9892 const size_t frame_size,
9893 const size_t hop_size,
9894 const bool pad_end =
true)
9905 "FFT::frame_signal");
9907 for (
size_t item = 0; item < offsets.
size(); ++item)
9909 const size_t offset = offsets[item];
9910 const size_t length = std::min(frame_size, signal.
size() -
offset);
9913 for (
size_t i = 0; i < frame_size; ++i)
9915 for (
size_t i = 0; i < length; ++i)
9916 frame(i) = signal[
offset + i];
9924 template <
typename Container>
9928 const size_t frame_size,
9929 const size_t hop_size,
9930 const bool pad_end =
true)
9943 const size_t hop_size,
9944 const bool pad_end =
true)
9955 const size_t frame_size,
9956 const size_t hop_size,
9957 const bool pad_end =
true)
9962 template <
typename SignalContainer,
typename WindowContainer>
9967 const size_t hop_size,
9968 const bool pad_end =
true)
9973 template <
typename Container>
9977 const size_t frame_size,
9978 const size_t hop_size,
9979 const bool pad_end =
true)
9999 const size_t chunk_size = 0)
10007 const size_t frame_size,
10017 const size_t frame_size,
10019 const size_t chunk_size = 0)
10024 template <
typename SignalContainer,
typename WindowContainer>
10034 template <
typename SignalContainer,
typename WindowContainer>
10041 const size_t chunk_size = 0)
10050 template <
typename Container>
10054 const size_t frame_size,
10060 template <
typename Container>
10065 const size_t frame_size,
10067 const size_t chunk_size = 0)
10083 const size_t hop_size,
10084 const size_t signal_length = 0)
10088 options.signal_length = signal_length;
10103 const size_t hop_size,
10104 const size_t signal_length = 0)
10112 const size_t frame_size,
10113 const size_t hop_size,
10114 const size_t signal_length = 0)
10126 const size_t hop_size,
10127 const size_t signal_length = 0,
10128 const size_t chunk_size = 0)
10132 options.signal_length = signal_length;
10146 const size_t hop_size,
10147 const size_t signal_length = 0,
10148 const size_t chunk_size = 0)
10163 const size_t frame_size,
10164 const size_t hop_size,
10165 const size_t signal_length = 0,
10166 const size_t chunk_size = 0)
10205 const size_t frame_size,
10219 const size_t chunk_size = 0)
10235 const size_t chunk_size = 0)
10244 const size_t frame_size,
10246 const size_t chunk_size = 0)
10284 > std::numeric_limits<size_t>::max() - block.size())
10285 <<
"FFT::STFTProcessor::process_block: pending buffer overflow";
10291 const size_t chunk_size,
10328 "FFT::STFTProcessor")),
10352 template <
typename WindowContainer>
10412 <<
"FFT::STFTProcessor::process_block: processor already flushed";
10414 if (block.is_empty())
10425 const size_t chunk_size = 0)
10429 <<
"FFT::STFTProcessor::process_block: processor already flushed";
10431 if (block.is_empty())
10495 template <
typename Container>
10503 template <
typename Container>
10508 const size_t chunk_size = 0)
10537 << ctx <<
": processor is not configured";
10573 > std::numeric_limits<size_t>::max()
10575 <<
"FFT::ISTFTProcessor::process_frame: frame offset overflow";
10579 <<
"FFT::ISTFTProcessor::process_frame: invalid overlap-add state";
10597 <<
"FFT::ISTFTProcessor::normalize_prefix: count " <<
count
10604 *
Real(256) * std::numeric_limits<Real>::epsilon();
10606 for (
size_t i = 0; i <
count; ++i)
10611 <<
"FFT::ISTFTProcessor: overlap-add normalization vanished at sample "
10642 size_t available = end - start;
10648 available = std::min(available,
10654 for (
size_t i = 0; i < available; ++i)
10670 <<
"FFT::ISTFTProcessor: invalid finalized-sample accounting";
10683 const size_t chunk_size)
10687 <<
"FFT::ISTFTProcessor::process_frame: processor already flushed";
10689 <<
"FFT::ISTFTProcessor::process_frame: frame size " <<
spectrum.
size()
10690 <<
" does not match configured FFT size " <<
fft_size_;
10727 "FFT::ISTFTProcessor");
10767 template <
typename AnalysisContainer,
typename SynthesisContainer>
10825 const size_t chunk_size = 0)
10842 const size_t chunk_size = 0)
10879 template <
typename Container>
10887 template <
typename Container>
10892 const size_t chunk_size = 0)
10907 << ctx <<
": batch size " <<
size
10908 <<
" does not match configured channel count " <<
processors_.size();
10917 <<
"FFT::BatchedSTFTProcessor: at least one channel is required";
10924 const size_t frame_size,
10929 template <
typename WindowContainer>
10948 <<
"FFT::BatchedSTFTProcessor::channel: index " << index
10949 <<
" out of range for " <<
processors_.size() <<
" channels";
10957 <<
"FFT::BatchedSTFTProcessor::channel: index " << index
10958 <<
" out of range for " <<
processors_.size() <<
" channels";
10982 const size_t chunk_size = 0)
10985 "FFT::BatchedSTFTProcessor::pprocess_block");
10991 [
this, &block, &pool, chunk_size, &
output](
const size_t i)
10994 processors_[i].pprocess_block(pool,
11020 [
this, &pool, chunk_size, &
output](
const size_t i)
11022 output(i) = processors_[i].pflush(pool, chunk_size);
11037 const size_t index)
11049 << ctx <<
": batch size " <<
size
11050 <<
" does not match configured channel count " <<
processors_.size();
11055 const size_t fft_size,
11068 const size_t fft_size,
11075 <<
"FFT::BatchedISTFTProcessor: at least one channel is required";
11078 <<
"FFT::BatchedISTFTProcessor: signal length count "
11092 const size_t fft_size,
11104 const size_t fft_size,
11117 const size_t fft_size,
11118 const size_t frame_size,
11129 const size_t fft_size,
11130 const size_t frame_size,
11141 template <
typename AnalysisContainer,
typename SynthesisContainer>
11145 const size_t fft_size,
11157 template <
typename AnalysisContainer,
typename SynthesisContainer>
11161 const size_t fft_size,
11183 <<
"FFT::BatchedISTFTProcessor::channel: index " << index
11184 <<
" out of range for " <<
processors_.size() <<
" channels";
11192 <<
"FFT::BatchedISTFTProcessor::channel: index " << index
11193 <<
" out of range for " <<
processors_.size() <<
" channels";
11207 "FFT::BatchedISTFTProcessor::process_block");
11217 const size_t chunk_size = 0)
11220 "FFT::BatchedISTFTProcessor::pprocess_block");
11225 [
this, &block, &pool, chunk_size, &
output](
const size_t i)
11228 processors_[i].pprocess_block(pool,
11252 [
this, &pool, chunk_size, &
output](
const size_t i)
11254 output(i) = processors_[i].pflush(pool, chunk_size);
11271 for (
size_t i = 0; i <
signals.size(); ++i)
11282 const size_t chunk_size = 0)
11297 for (
size_t i = 0; i <
signals.size(); ++i)
11306 const size_t frame_size,
11316 const size_t frame_size,
11318 const size_t chunk_size = 0)
11340 <<
"FFT::batched_istft: signal_lengths size "
11342 <<
" does not match batch size " <<
spectrograms.size();
11368 const size_t chunk_size = 0)
11375 <<
"FFT::pbatched_istft: signal_lengths size "
11377 <<
" does not match batch size " <<
spectrograms.size();
11436 const size_t chunk_size = 0)
11450 const size_t frame_size,
11462 const size_t frame_size,
11465 const size_t chunk_size = 0)
11488 << ctx <<
": filter is not configured";
11496 for (
size_t i = 0; i <
order; ++i)
11524 :
LFilter(coeffs.numerator, coeffs.denominator)
11540 template <
typename NumContainer,
typename DenContainer>
11584 <<
"FFT::LFilter::set_state: state size " <<
new_state.size()
11585 <<
" does not match filter order " <<
state_.
size();
11598 "FFT::LFilter::filter",
11603 template <
typename Container>
11624 << ctx <<
": filter is not configured";
11632 for (
size_t i = 0; i <
coeffs_.size(); ++i)
11635 for (
size_t j = 0; j <
state.
size(); ++j)
11652 <<
"FFT::SOSFilter: at least one biquad section is required";
11655 for (
size_t i = 0; i <
sections_.size(); ++i)
11658 "FFT::SOSFilter"));
11665 template <
typename SectionsContainer>
11688 <<
"FFT::SOSFilter::state: index " << index
11689 <<
" out of range for " <<
states_.
size() <<
" sections";
11697 for (
size_t i = 0; i <
states_.size(); ++i)
11698 for (
size_t j = 0; j <
states_[i].size(); ++j)
11708 for (
size_t i = 0; i <
coeffs_.size(); ++i)
11715 "FFT::SOSFilter::filter",
11723 template <
typename Container>
11742 << ctx <<
": batch size " <<
count
11743 <<
" does not match configured channel count " <<
filters_.size();
11752 <<
"FFT::LFilterBank: at least one channel is required";
11768 template <
typename NumContainer,
typename DenContainer>
11787 <<
"FFT::LFilterBank::channel: index " << index
11788 <<
" out of range for " <<
filters_.size() <<
" channels";
11796 <<
"FFT::LFilterBank::channel: index " << index
11797 <<
" out of range for " <<
filters_.size() <<
" channels";
11803 for (
size_t i = 0; i <
filters_.size(); ++i)
11814 template <
typename Container>
11828 for (
size_t i = 0; i <
filters_.size(); ++i)
11836 const size_t chunk_size = 0)
11845 output(i) = filters_[i].filter(signals[i]);
11862 << ctx <<
": batch size " <<
count
11863 <<
" does not match configured channel count " <<
filters_.size();
11871 <<
"FFT::SOSFilterBank: at least one channel is required";
11877 template <
typename SectionsContainer>
11893 <<
"FFT::SOSFilterBank::channel: index " << index
11894 <<
" out of range for " <<
filters_.size() <<
" channels";
11902 <<
"FFT::SOSFilterBank::channel: index " << index
11903 <<
" out of range for " <<
filters_.size() <<
" channels";
11909 for (
size_t i = 0; i <
filters_.size(); ++i)
11920 template <
typename Container>
11934 for (
size_t i = 0; i <
filters_.size(); ++i)
11942 const size_t chunk_size = 0)
11951 output(i) = filters_[i].filter(signals[i]);
11965 const IIRCoefficients coeffs =
11969 coeffs.denominator,
11990 template <
typename SignalContainer,
typename NumContainer,
typename DenContainer>
12016 const IIRCoefficients coeffs =
12020 <<
"FFT::batched_lfilter: initial state count "
12022 <<
" does not match channel count " <<
signals.size();
12025 for (
size_t i = 0; i <
signals.size(); ++i)
12028 coeffs.denominator,
12030 "FFT::batched_lfilter");
12040 const size_t chunk_size = 0)
12045 const IIRCoefficients coeffs =
12049 <<
"FFT::pbatched_lfilter: initial state count "
12051 <<
" does not match channel count " <<
signals.size();
12059 output(i) = iir_filter_impl(signals[i],
12061 coeffs.denominator,
12062 initial_states.is_empty() ? Array<Real>() : initial_states[i],
12063 "FFT::pbatched_lfilter");
12069 template <
typename SignalsContainer,
typename NumContainer,
typename DenContainer>
12080 for (
const auto & signal:
signals)
12096 template <
typename SignalContainer,
typename SectionsContainer>
12115 for (
size_t i = 0; i <
signals.size(); ++i)
12124 const size_t chunk_size = 0)
12135 output(i) = sosfilt(signals[i], sections);
12141 template <
typename SignalsContainer,
typename SectionsContainer>
12149 for (
const auto & signal:
signals)
12159 const bool whole =
false)
12168 const bool whole =
false)
12176 const bool whole =
false)
12184 const bool whole =
false)
12192 const bool whole =
false)
12195 <<
"FFT::freqz: at least one biquad section is required";
12197 <<
"FFT::freqz: number of frequency samples must be positive";
12210 for (
size_t j = 0; j <
sections.size(); ++j)
12214 "FFT::freqz(SOS)");
12215 output.omega(i) = omega;
12216 output.response(i) = response;
12222 template <
typename NumContainer,
typename DenContainer>
12228 const bool whole =
false)
12236 template <
typename SectionsContainer>
12241 const bool whole =
false)
12300 for (
size_t i = 0; i <
sections.size(); ++i)
12309 for (
size_t i = 0; i <
sections.size(); ++i)
12314 template <
typename Container>
12322 template <
typename Container>
12330 template <
typename SectionsContainer>
12338 template <
typename SectionsContainer>
12379 template <
typename NumContainer,
typename DenContainer>
12389 template <
typename SectionsContainer>
12430 template <
typename NumContainer,
typename DenContainer>
12440 template <
typename SectionsContainer>
12473 Real margin = std::numeric_limits<Real>::infinity();
12474 for (
size_t i = 0; i <
sections.size(); ++i)
12479 template <
typename Container>
12487 template <
typename SectionsContainer>
12499 const Real tol =
Real(1024) * std::numeric_limits<Real>::epsilon();
12518 for (
size_t i = 0; i <
sections.size(); ++i)
12524 template <
typename Container>
12532 template <
typename SectionsContainer>
12544 const Real tolerance)
12547 for (
size_t i = 0; i < pairs.
size(); ++i)
12548 if (pairs[i].is_cancellation(tolerance))
12556 const Real tolerance)
12559 poles(denominator),
12565 const Real tolerance)
12574 const Real tolerance)
12583 const Real tolerance)
12590 template <
typename NumContainer,
typename DenContainer>
12595 const Real tolerance)
12602 template <
typename SectionsContainer>
12606 const Real tolerance)
12615 const Real tol =
Real(1024) * std::numeric_limits<Real>::epsilon();
12619 <<
"FFT::validate_stable: stability margin " <<
margin
12620 <<
" is not strictly positive";
12621 for (
size_t i = 0; i < roots.
size(); ++i)
12623 <<
"FFT::validate_stable: unstable pole at " << roots[i]
12624 <<
" with magnitude " << std::abs(roots[i]);
12642 for (
size_t i = 0; i <
sections.size(); ++i)
12646 template <
typename Container>
12654 template <
typename SectionsContainer>
12668 <<
"FFT::validate_stable: stability margin " <<
margin
12669 <<
" is smaller than required margin " <<
min_margin;
12693 <<
"FFT::validate_stable: SOS stability margin "
12695 <<
" is smaller than required margin " <<
min_margin;
12699 template <
typename Container>
12708 template <
typename SectionsContainer>
12720 const Real tolerance)
12723 for (
size_t i = 0; i < pairs.
size(); ++i)
12724 if (pairs[i].is_cancellation(tolerance))
12726 <<
"FFT::validate_no_near_pole_zero_cancellation: zero "
12727 << pairs[i].zero <<
" and pole " << pairs[i].pole
12728 <<
" are only " << pairs[i].distance <<
" apart";
12734 const Real tolerance)
12737 poles(denominator),
12743 const Real tolerance)
12752 const Real tolerance)
12761 const Real tolerance)
12768 template <
typename NumContainer,
typename DenContainer>
12773 const Real tolerance)
12780 template <
typename SectionsContainer>
12784 const Real tolerance)
12808 const bool whole =
false)
12814 "FFT::group_delay");
12821 const bool whole =
false)
12827 "FFT::phase_delay");
12833 const bool whole =
false)
12839 "FFT::group_delay");
12845 const bool whole =
false)
12851 "FFT::phase_delay");
12857 const bool whole =
false)
12865 const bool whole =
false)
12873 const bool whole =
false)
12884 const bool whole =
false)
12895 const bool whole =
false)
12900 "FFT::group_delay(SOS)");
12906 const bool whole =
false)
12913 "FFT::phase_delay(SOS)");
12916 *
Real(128) * std::numeric_limits<Real>::epsilon();
12918 for (
size_t i = 0; i < response.
omega.
size(); ++i)
12924 template <
typename NumContainer,
typename DenContainer>
12930 const bool whole =
false)
12938 template <
typename NumContainer,
typename DenContainer>
12944 const bool whole =
false)
12952 template <
typename Container>
12957 const bool whole =
false)
12962 template <
typename Container>
12967 const bool whole =
false)
12972 template <
typename SectionsContainer>
12977 const bool whole =
false)
12982 template <
typename SectionsContainer>
12987 const bool whole =
false)
13010 const bool whole =
false)
13018 [&](
const Real omega)
13023 "FFT::phase_margin");
13031 const bool whole =
false)
13039 [&](
const Real omega)
13044 "FFT::gain_margin");
13051 const bool whole =
false)
13059 const bool whole =
false)
13067 const bool whole =
false)
13078 const bool whole =
false)
13089 const bool whole =
false)
13094 [&](
const Real omega)
13098 "FFT::phase_margin(SOS)");
13105 const bool whole =
false)
13110 [&](
const Real omega)
13114 "FFT::gain_margin(SOS)");
13118 template <
typename NumContainer,
typename DenContainer>
13124 const bool whole =
false)
13132 template <
typename NumContainer,
typename DenContainer>
13138 const bool whole =
false)
13146 template <
typename SectionsContainer>
13151 const bool whole =
false)
13156 template <
typename SectionsContainer>
13161 const bool whole =
false)
13175 "FFT::bilinear_transform");
13178 template <
typename NumContainer,
typename DenContainer>
13197 "FFT::butterworth_lowpass"),
13202 "FFT::butterworth_lowpass");
13212 "FFT::butterworth_highpass"),
13217 "FFT::butterworth_highpass");
13231 "FFT::butterworth_bandpass");
13235 "FFT::butterworth_bandpass"),
13236 {
Real(1),
Real(0), center * center},
13239 "FFT::butterworth_bandpass");
13254 "FFT::butterworth_bandstop");
13266 "FFT::chebyshev1_lowpass"),
13271 "FFT::chebyshev1_lowpass");
13283 "FFT::chebyshev1_highpass"),
13288 "FFT::chebyshev1_highpass");
13303 "FFT::chebyshev1_bandpass");
13308 "FFT::chebyshev1_bandpass"),
13309 {
Real(1),
Real(0), center * center},
13312 "FFT::chebyshev1_bandpass");
13326 "FFT::chebyshev1_bandstop"),
13331 "FFT::chebyshev1_bandstop");
13337 const Real attenuation_db,
13343 "FFT::chebyshev2_lowpass"),
13347 "FFT::chebyshev2_lowpass");
13353 const Real attenuation_db,
13359 "FFT::chebyshev2_highpass"),
13363 "FFT::chebyshev2_highpass");
13369 const Real attenuation_db,
13378 "FFT::chebyshev2_bandpass");
13383 "FFT::chebyshev2_bandpass"),
13384 {
Real(1),
Real(0), center * center},
13387 "FFT::chebyshev2_bandpass");
13393 const Real attenuation_db,
13402 "FFT::chebyshev2_bandstop");
13407 "FFT::chebyshev2_bandstop"),
13409 {
Real(1),
Real(0), center * center},
13411 "FFT::chebyshev2_bandstop");
13425 "FFT::bessel_lowpass"),
13429 "FFT::bessel_lowpass");
13439 "FFT::bessel_highpass"),
13443 "FFT::bessel_highpass");
13457 "FFT::bessel_bandpass");
13461 "FFT::bessel_bandpass"),
13462 {
Real(1),
Real(0), center * center},
13465 "FFT::bessel_bandpass");
13480 "FFT::bessel_bandstop");
13487 const Real attenuation_db,
13494 "FFT::elliptic_lowpass"),
13498 "FFT::elliptic_lowpass");
13505 const Real attenuation_db,
13512 "FFT::elliptic_highpass"),
13516 "FFT::elliptic_highpass");
13523 const Real attenuation_db,
13532 "FFT::elliptic_bandpass");
13538 "FFT::elliptic_bandpass"),
13539 {
Real(1),
Real(0), center * center},
13542 "FFT::elliptic_bandpass");
13549 const Real attenuation_db,
13558 "FFT::elliptic_bandstop");
13564 "FFT::elliptic_bandstop"),
13566 {
Real(1),
Real(0), center * center},
13568 "FFT::elliptic_bandstop");
13580 const size_t block_size = 0)
13585 if (coeffs.
size() == 1)
13588 const Real gain = coeffs[0] * coeffs[0];
13589 for (
size_t i = 0; i < signal.
size(); ++i)
13590 output(i) = signal[i] * gain;
13611 const size_t block_size = 0,
13612 const size_t chunk_size = 0)
13617 if (coeffs.
size() == 1)
13620 const Real gain = coeffs[0] * coeffs[0];
13621 for (
size_t i = 0; i < signal.
size(); ++i)
13622 output(i) = signal[i] * gain;
13640 template <
typename SignalContainer,
typename CoeffContainer>
13645 const size_t block_size = 0)
13650 template <
typename SignalContainer,
typename CoeffContainer>
13656 const size_t block_size = 0,
13657 const size_t chunk_size = 0)
13681 "FFT::filtfilt(IIR)");
13708 template <
typename SignalContainer,
typename NumContainer,
typename DenContainer>
13722 template <
typename SignalContainer>
13731 template <
typename SignalContainer>
13740 template <
typename SignalContainer,
typename SectionsContainer>
13771 <<
"FFT::OverlapAdd: kernel size must be positive";
13781 const size_t length)
const
13787 for (
size_t i = 0; i < length; ++i)
13799 const size_t chunk_size)
const
13832 const size_t chunk_size)
13834 if (block.is_empty())
13838 <<
"FFT::OverlapAdd::process_block: block size " << block.size()
13839 <<
" exceeds configured block size " <<
block_size_;
13843 if (pool !=
nullptr)
13855 for (
size_t i = 0; i < block.size(); ++i)
13866 const size_t src = block.size() + i;
13882 const size_t chunk_size)
const
13888 > std::numeric_limits<size_t>::max()
13890 <<
"FFT::OverlapAdd::convolve: product size exceeds size_t capacity";
13893 for (
size_t i = 0; i <
output.size(); ++i)
13901 if (pool !=
nullptr)
13912 const size_t valid = length +
kernel_.
size() - 1;
13913 for (
size_t i = 0; i < valid; ++i)
13926 <<
"FFT::OverlapAdd: kernel must be non-empty";
13930 <<
"FFT::OverlapAdd: block size must be positive";
13933 > std::numeric_limits<size_t>::max()
13935 <<
"FFT::OverlapAdd: FFT size overflow";
13970 const size_t chunk_size = 0)
const
13978 if (block.is_empty())
13988 for (
size_t i = 0; i < length; ++i)
13992 for (
size_t i = 0; i < partial.
size(); ++i)
14001 const size_t chunk_size = 0)
14003 if (block.is_empty())
14013 for (
size_t i = 0; i < length; ++i)
14017 for (
size_t i = 0; i < partial.
size(); ++i)
14038 template <
typename Container>
14046 template <
typename Container>
14050 const size_t chunk_size = 0)
const
14055 template <
typename Container>
14063 template <
typename Container>
14067 const size_t chunk_size = 0)
14096 <<
"FFT::OverlapAddBank: kernel size must be positive";
14105 << ctx <<
": batch size " <<
count
14106 <<
" does not match configured channel count " <<
overlaps_.size();
14114 for (
size_t i = 0; i <
signals.size(); ++i)
14123 const size_t length)
14127 for (
size_t i = 0; i < length; ++i)
14135 const size_t offset = 0,
14136 const size_t length = std::numeric_limits<size_t>::max())
const
14139 <<
"FFT::OverlapAddBank::build_signal_block: offset " <<
offset
14140 <<
" exceeds signal size " << signal.
size();
14144 <<
"FFT::OverlapAddBank::build_signal_block: block size "
14160 for (
size_t channel = 0; channel < block.size(); ++channel)
14170 const size_t offset = 0)
const
14173 for (
size_t item = 0; item <
channels.size(); ++item)
14186 const size_t chunk_size)
const
14197 for (
size_t channel = 0; channel <
spectra.size(); ++channel)
14205 for (
size_t i = 0; i <
overlaps_[channel].size(); ++i)
14213 for (
size_t channel = 0; channel <
overlaps_.size(); ++channel)
14220 const size_t chunk_size)
14229 for (
size_t channel = 0; channel < block.size(); ++channel)
14230 if (
not block[channel].is_empty())
14269 for (
size_t i = 0; i <
overlaps_[channel].size(); ++i)
14288 const size_t chunk_size)
const
14294 for (
size_t channel = 0; channel <
signals.size(); ++channel)
14296 if (
signals[channel].is_empty())
14300 for (
size_t i = 0; i <
output[channel].size(); ++i)
14311 for (
size_t channel = 0; channel <
signals.size(); ++channel)
14316 const size_t length =
14348 for (
size_t i = 0; i < valid; ++i)
14363 <<
"FFT::OverlapAddBank: at least one channel is required";
14365 <<
"FFT::OverlapAddBank: kernel must be non-empty";
14369 <<
"FFT::OverlapAddBank: block size must be positive";
14371 > std::numeric_limits<size_t>::max()
14373 <<
"FFT::OverlapAddBank: FFT size overflow";
14380 for (
size_t channel = 0; channel <
num_channels; ++channel)
14422 const size_t chunk_size = 0)
const
14434 for (
size_t channel = 0; channel <
emitted.size(); ++channel)
14440 for (
size_t channel = 0; channel < block.size(); ++channel)
14444 const size_t length =
14455 for (
size_t channel = 0; channel <
emitted.size(); ++channel)
14456 for (
size_t i = 0; i < partial[channel].
size(); ++i)
14457 emitted(channel).append(partial[channel][i]);
14466 const size_t chunk_size = 0)
14471 for (
size_t channel = 0; channel <
emitted.size(); ++channel)
14477 for (
size_t channel = 0; channel < block.size(); ++channel)
14481 const size_t length =
14492 for (
size_t channel = 0; channel <
emitted.size(); ++channel)
14493 for (
size_t i = 0; i < partial[channel].
size(); ++i)
14494 emitted(channel).append(partial[channel][i]);
14504 for (
size_t channel = 0; channel <
overlaps_.size(); ++channel)
14517 auto flush_one = [
this, &tail](
const size_t channel)
14527 for (
size_t channel = 0; channel <
overlaps_.size(); ++channel)
14558 <<
"FFT::OverlapSave: kernel size must be positive";
14605 <<
"FFT::OverlapSave::process_block: block size " << chunk.
size()
14606 <<
" exceeds configured block size " <<
block_size_;
14671 <<
"FFT::OverlapSave: kernel must be non-empty";
14675 <<
"FFT::OverlapSave: block size must be positive";
14677 > std::numeric_limits<size_t>::max()
14679 <<
"FFT::OverlapSave: FFT size overflow";
14727 if (block.is_empty())
14797 template <
typename Container>
14806 template <
typename Container>
14841 <<
"FFT::PartitionedConvolver: kernel size must be positive";
14895 for (
size_t i = 0; i <
padded.size(); ++i)
14928 <<
"FFT::PartitionedConvolver: kernel must be non-empty";
14933 <<
"FFT::PartitionedConvolver: partition size must be positive";
14941 const size_t num_partitions =
14949 const size_t length =
14951 for (
size_t i = 0; i < length; ++i)
14984 if (block.is_empty())
15045 template <
typename Container>
15053 template <
typename Container>
15066 const size_t block_size = 0)
15079 const size_t block_size = 0,
15080 const size_t chunk_size = 0)
15091 const size_t block_size = 0)
15103 const size_t block_size = 0,
15104 const size_t chunk_size = 0)
15116 const size_t block_size = 0)
15127 const size_t partition_size = 0)
15135 template <
typename Container1,
typename Container2>
15140 const size_t block_size = 0)
15148 template <
typename Container1,
typename Container2>
15154 const size_t block_size = 0,
15155 const size_t chunk_size = 0)
15165 template <
typename Container1,
typename Container2>
15170 const size_t block_size = 0)
15178 template <
typename Container1,
typename Container2>
15183 const size_t partition_size = 0)
15218 const size_t chunk_size = 0)
15224 template <
typename Container1,
typename Container2>
15233 template <
typename Container1,
typename Container2>
15237 const size_t chunk_size = 0)
15264 const size_t chunk_size = 0)
15270 template <
typename Container1,
typename Container2>
15279 template <
typename Container1,
typename Container2>
15283 const size_t chunk_size = 0)
Exception handling system with formatted messages for Aleph-w.
#define ah_out_of_range_error_if(C)
Throws std::out_of_range if condition holds.
#define ah_domain_error_if(C)
Throws std::domain_error if condition holds.
#define ah_runtime_error_if(C)
Throws std::runtime_error if condition holds.
#define ah_invalid_argument_if(C)
Throws std::invalid_argument if condition holds.
size_t size_t int32_t value
Simple dynamic array with automatic resizing and functional operations.
static Array create(size_t n)
Create an array with n logical elements.
constexpr size_t size() const noexcept
Return the number of elements stored in the stack.
void empty() noexcept
Empties the container.
constexpr bool is_empty() const noexcept
Checks if the container is empty.
T & append(const T &data)
Append a copy of data
void reserve(size_t cap)
Reserves cap cells into the array.
Stateful multichannel ISTFT processor with one synthesizer per channel.
ISTFTProcessor & channel(const size_t index)
void validate_batch_size(const size_t size, const char *ctx) const
BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const size_t frame_size, const ISTFTOptions &options)
BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const Array< Real > &window, const ISTFTOptions &options)
size_t num_channels() const noexcept
static ISTFTOptions channel_options(const ISTFTOptions &options, const Array< size_t > &signal_lengths, const size_t index)
BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const Array< Real > &window, const ISTFTOptions &options, const Array< size_t > &signal_lengths)
and Is_Real_Container< SynthesisContainer > BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const AnalysisContainer &analysis_window, const SynthesisContainer &synthesis_window, const ISTFTOptions &options, const Array< size_t > &signal_lengths)
Array< Array< Real > > pprocess_block(ThreadPool &pool, const Array< Array< Array< Complex > > > &block, const size_t chunk_size=0)
const ISTFTProcessor & channel(const size_t index) const
BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const size_t frame_size, const ISTFTOptions &options, const Array< size_t > &signal_lengths)
Array< Array< Real > > flush()
and Is_Real_Container< SynthesisContainer > BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const AnalysisContainer &analysis_window, const SynthesisContainer &synthesis_window, const ISTFTOptions &options)
BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options)
Array< Array< Real > > pflush(ThreadPool &pool, const size_t chunk_size=0)
Array< ISTFTProcessor > processors_
Array< Array< Real > > process_block(const Array< Array< Array< Complex > > > &block)
BatchedISTFTProcessor(const size_t num_channels, const size_t fft_size, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options, const Array< size_t > &signal_lengths)
Stateful multichannel STFT processor with one analyzer per channel.
Array< Array< Array< Complex > > > pflush(ThreadPool &pool, const size_t chunk_size=0)
BatchedSTFTProcessor(const size_t num_channels, const size_t frame_size, const STFTOptions &options)
Array< STFTProcessor > processors_
Array< Array< Array< Complex > > > process_block(const Array< Array< Real > > &block)
const STFTProcessor & channel(const size_t index) const
size_t num_channels() const noexcept
Array< Array< Array< Complex > > > flush()
STFTProcessor & channel(const size_t index)
void validate_batch_size(const size_t size, const char *ctx) const
Array< Array< Array< Complex > > > pprocess_block(ThreadPool &pool, const Array< Array< Real > > &block, const size_t chunk_size=0)
BatchedSTFTProcessor(const size_t num_channels, const Array< Real > &window, const STFTOptions &options)
BatchedSTFTProcessor(const size_t num_channels, const WindowContainer &window, const STFTOptions &options)
Stateful ISTFT processor for chunked frame-by-frame synthesis.
Array< Real > analysis_window_
ISTFTProcessor(const size_t fft_size, const size_t frame_size, const ISTFTOptions &options)
Construct an inverse STFT processor with a Hann window.
Array< Real > synthesis_window_
bool centered() const noexcept
size_t left_trim_remaining_
Array< Real > pprocess_frame(ThreadPool &pool, const Container &spectrum, const size_t chunk_size=0)
Array< Real > pprocess_block(ThreadPool &pool, const Array< Array< Complex > > &spectrogram_block, const size_t chunk_size=0)
size_t fft_size() const noexcept
Array< Real > drain_ready_samples(const bool final_flush)
size_t hop_size() const noexcept
Array< Real > process_frame_impl(const Array< Complex > &spectrum, ThreadPool *pool, const size_t chunk_size)
Array< Real > process_frame(const Array< Complex > &spectrum)
and Is_Real_Container< SynthesisContainer > ISTFTProcessor(const size_t fft_size, const AnalysisContainer &analysis_window, const SynthesisContainer &synthesis_window, const ISTFTOptions &options)
Construct an inverse STFT processor from any real containers.
Array< Real > emit_samples(const Array< Real > &normalized, const bool final_flush)
Array< Real > pprocess_frame(ThreadPool &pool, const Array< Complex > &spectrum, const size_t chunk_size=0)
void require_configured(const char *ctx) const
Array< Real > pflush(ThreadPool &, const size_t=0)
Array< Real > process_block(const Array< Array< Complex > > &spectrogram_block)
Array< Real > pending_norm_
Array< Real > pending_output_
void accumulate_frame(const Array< Real > &frame)
ISTFTProcessor(const size_t fft_size, const Array< Real > &window, const ISTFTOptions &options)
Construct an inverse STFT processor with a single window.
Array< Real > process_frame(const Container &spectrum)
ISTFTProcessor(const size_t fft_size, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options)
Construct an inverse STFT processor with custom windows.
size_t emitted_raw_samples_
ISTFTProcessor()=default
Default constructor for an unconfigured inverse processor.
void ensure_pending_size(const size_t size)
bool finalized() const noexcept
Array< Real > normalize_prefix(const size_t count)
size_t frame_size() const noexcept
size_t emitted_signal_samples_
bool configured() const noexcept
Stateful bank of direct-form II transposed IIR filters.
LFilterBank(const size_t num_channels, const IIRCoefficients &coeffs)
void validate_channel_count(const size_t count, const char *ctx) const
Rationale: Internal safety guard.
LFilter & channel(const size_t index)
size_t num_channels() const noexcept
Array< Real > filter_channel(const size_t index, const Array< Real > &signal)
Array< LFilter > filters_
and Is_Real_Container< DenContainer > LFilterBank(const size_t num_channels, const NumContainer &numerator, const DenContainer &denominator)
LFilterBank(const size_t num_channels, const Array< Real > &numerator, const Array< Real > &denominator)
Array< Array< Real > > pfilter(ThreadPool &pool, const Array< Array< Real > > &signals, const size_t chunk_size=0)
Array< Array< Real > > filter(const Array< Array< Real > > &signals)
LFilterBank(const size_t num_channels, const BiquadSection §ion)
const LFilter & channel(const size_t index) const
Array< Real > filter_channel(const size_t index, const Container &signal)
Stateful direct-form II transposed IIR filter.
LFilter(const IIRCoefficients &coeffs)
Construct a linear filter from IIRCoefficients.
bool configured() const noexcept
Returns whether the filter has been initialized with coefficients.
and Is_Real_Container< DenContainer > LFilter(const NumContainer &numerator, const DenContainer &denominator)
Construct a linear filter from any real containers.
void require_configured(const char *ctx) const
const Array< Real > & state() const noexcept
Returns the current internal delay-line state.
const IIRCoefficients & coefficients() const noexcept
Returns the active normalized filter coefficients.
void reset()
Resets the delay-line state to zero.
LFilter(const Array< Real > &numerator, const Array< Real > &denominator)
Construct a linear filter from numerator and denominator arrays.
Array< Real > filter(const Container &signal)
Processes a real-valued container through the stateful filter.
void set_state(const Array< Real > &new_state)
Manually overrides the delay-line state.
LFilter(const BiquadSection §ion)
Construct a linear filter from a single BiquadSection.
size_t order() const noexcept
Returns the filter order (number of feedback coefficients - 1).
LFilter()=default
Default constructor for an unconfigured LFilter.
Array< Real > filter(const Array< Real > &signal)
Processes a signal block through the stateful filter.
Multichannel overlap-add convolver with one shared kernel FFT.
Array< bool > has_pending_tail_
void validate_channel_count(const size_t count, const char *ctx) const
Rationale: Internal safety guard.
Array< Array< Real > > convolve(const Array< Array< Real > > &signals) const
const Array< Real > & kernel() const noexcept
Array< Array< Real > > convolve_impl(const Array< Array< Real > > &signals, ThreadPool *pool, const size_t chunk_size) const
static Array< Real > slice_chunk(const Array< Real > &input, const size_t offset, const size_t length)
Returns a contiguous slice of input.
Array< Array< Real > > pflush(ThreadPool &pool, const size_t chunk_size=0)
Array< Array< Real > > process_chunk_batch_impl(const Array< Array< Real > > &block, ThreadPool *pool, const size_t chunk_size)
Array< Array< Real > > process_block(const Array< Array< Real > > &block)
Array< Array< Real > > overlaps_
static size_t max_batch_length(const Array< Array< Real > > &signals) noexcept
Returns the maximum signal length in a batch.
size_t block_size() const noexcept
size_t overlap_size() const noexcept
Array< Array< Real > > pconvolve(ThreadPool &pool, const Array< Array< Real > > &signals, const size_t chunk_size=0) const
Array< Complex > kernel_spectrum_
void pointwise_multiply_batch(Array< Array< Complex > > &spectra, ThreadPool *pool, const size_t chunk_size) const
Rationale: Parallel pointwise multiplication of multiple signal spectra with the shared kernel spectr...
Array< Array< Real > > pprocess_block(ThreadPool &pool, const Array< Array< Real > > &block, const size_t chunk_size=0)
void clear_channel_overlap(const size_t channel)
Rationale: Resets the overlap buffer for a single channel.
static size_t default_block_size(const size_t kernel_size)
Rationale: Selects a block size that is at least as large as the kernel to ensure efficient FFT proce...
Array< Array< Complex > > build_signal_batch(const Array< Array< Real > > &source, const Array< size_t > &channels, const Array< size_t > &lengths, const size_t offset=0) const
implementation helper for preparing a subset of signal blocks.
OverlapAddBank(const size_t num_channels, const Array< Real > &kernel, const size_t block_size=0)
size_t fft_size() const noexcept
size_t num_channels() const noexcept
Array< Complex > build_signal_block(const Array< Real > &signal, const size_t offset=0, const size_t length=std::numeric_limits< size_t >::max()) const
implementation helper for preparing a signal block for FFT.
Array< Array< Complex > > build_signal_batch(const Array< Array< Real > > &block) const
implementation helper for preparing a batch of signal blocks.
Array< Array< Real > > flush()
Reusable overlap-add convolver for long real sequences.
Array< Real > convolve(const Container &signal) const
void clear_overlap()
Rationale: Resets the overlap buffer between independent convolution runs.
Array< Real > convolve_impl(const Array< Real > &signal, ThreadPool *pool, const size_t chunk_size) const
implementation of one-shot linear convolution.
Array< Complex > kernel_spectrum_
Array< Real > process_block(const Array< Real > &block)
size_t block_size() const noexcept
const Array< Real > & kernel() const noexcept
Array< Real > pprocess_block(ThreadPool &pool, const Array< Real > &block, const size_t chunk_size=0)
Array< Real > process_chunk_impl(const Array< Real > &block, ThreadPool *pool, const size_t chunk_size)
Rationale: Implements the core Overlap-Add logic for a single block:
Array< Real > pprocess_block(ThreadPool &pool, const Container &block, const size_t chunk_size=0)
OverlapAdd(const Array< Real > &kernel, const size_t block_size=0)
Array< Real > pconvolve(ThreadPool &pool, const Container &signal, const size_t chunk_size=0) const
Array< Real > convolve(const Array< Real > &signal) const
size_t overlap_size() const noexcept
Array< Real > process_block(const Container &block)
size_t fft_size() const noexcept
static size_t default_block_size(const size_t kernel_size)
Array< Real > pconvolve(ThreadPool &pool, const Array< Real > &signal, const size_t chunk_size=0) const
void pointwise_multiply(Array< Complex > &spectrum, ThreadPool *pool, const size_t chunk_size) const
Rationale: Applies the precomputed FIR kernel in the frequency domain via element-wise multiplication...
Array< Complex > build_signal_block(const Array< Real > &signal, const size_t offset, const size_t length) const
Rationale: Prepares a chunk of signal for FFT by zero-padding it to the internal fft_size().
Reusable overlap-save convolver for streaming real FIR filtering.
Array< Real > process_block(const Array< Real > &block)
Processes an input block and returns the corresponding filtered output.
void update_history(const Array< Real > &padded_chunk)
Rationale: Maintains the sliding window of historical samples required to correctly compute the linea...
size_t fft_size() const noexcept
Returns the FFT size used internally (block_size + kernel_size - 1 padded to power of 2).
Array< Real > process_block(const Container &block)
Processes an input block from a real-valued container.
void clear_history()
Rationale: Resets the state for a new filtering process.
Array< Complex > kernel_spectrum_
void reset()
Resets the internal filter state and history.
static size_t default_block_size(const size_t kernel_size)
Rationale: Heuristic for choosing a block size that provides a good balance between algorithmic laten...
Array< Real > process_chunk_impl(const Array< Real > &chunk)
Rationale: Implements the core Overlap-Save logic for one block:
Array< Complex > build_segment_spectrum(const Array< Real > &chunk) const
Rationale: Prepares a frequency-domain segment by combining historical samples with the current signa...
Array< Real > flush()
Flushes all buffered samples and returns the final tail of the convolution.
const Array< Real > & kernel() const noexcept
Returns the current FIR kernel.
size_t overlap_size() const noexcept
Returns the internal overlap history size.
size_t block_size() const noexcept
Returns the processing block size (the hop size).
Array< Real > convolve(const Container &signal)
Convolves an input signal from a real-valued container.
size_t emitted_output_size_
OverlapSave(const Array< Real > &kernel, const size_t block_size=0)
Array< Real > convolve(const Array< Real > &signal)
Convolves the input signal with the kernel (one-shot convenience).
Array< Real > padded_chunk_copy(const Array< Real > &chunk) const
Internal helper for zero-padding chunks.
Uniform partitioned FIR convolver for low-latency streaming.
Array< Real > convolve(const Container &signal)
Array< Real > process_block(const Container &block)
size_t fft_size() const noexcept
const Array< Real > & kernel() const noexcept
void clear_state()
Rationale: Resets all delay-line history and overlap buffers.
PartitionedConvolver(const Array< Real > &kernel, const size_t partition_size=0)
Array< Array< Complex > > input_history_
Array< Real > padded_partition(const Array< Real > &chunk) const
implementation helper for preparing a padded partition.
Array< Real > process_partition_impl(const Array< Real > &chunk)
Rationale: Implements the core frequency-domain partitioned convolution logic:
Array< Real > convolve(const Array< Real > &signal)
Array< Real > process_block(const Array< Real > &block)
size_t emitted_output_size_
static size_t default_partition_size(const size_t kernel_size)
Rationale: Selects a partition size that minimizes latency while maintaining FFT efficiency (capped a...
size_t partition_size() const noexcept
Array< Complex > zero_spectrum() const
implementation helper for creating zero-filled complex spectra.
Array< Array< Complex > > kernel_partitions_
Precomputed FFT plan for repeated transforms of the same size.
Array< Complex > rfft(const Array< Real > &input) const
Computes the Real-to-Complex FFT (RFFT).
std::shared_ptr< const Plan > half_plan_
void initialize_mixed_radix_plan()
void apply_bit_reversal(Array< Complex > &a) const noexcept
Array< Array< Complex > > ptransformed_batch(ThreadPool &pool, const Array< Array< Complex > > &input, const bool invert=false, const size_t chunk_size=0, const bool prefer_simd=true) const
Returns the parallel batch FFT/IFFT as a new array of arrays.
Array< Complex > ptransformed(ThreadPool &pool, const Array< Complex > &input, const bool invert=false, const size_t chunk_size=0) const
Parallel transform returning a new array.
void transform_batch(Array< Array< Complex > > &batch, const bool invert, const bool prefer_simd=true) const
In-place batch transform for equal-length complex inputs.
Array< Real > inverse_transform_real(const Array< Complex > &input) const
Returns the IFFT projected back to real values.
void initialize_power_of_two_plan()
void initialize_bluestein_plan()
Plan(const size_t n)
Constructs a plan for transforms of size n.
Array< Complex > prfft(ThreadPool &pool, const Array< Real > &input, const size_t chunk_size=0) const
Parallel Real-to-Complex FFT (RFFT).
Array< Real > pirfft(ThreadPool &pool, const Array< Complex > &spectrum, const size_t chunk_size=0) const
Parallel Complex-to-Real Inverse FFT (IRFFT).
Complex root_for_length(const size_t length, const size_t exponent, const bool invert) const
void transform(Array< Complex > &a, const bool invert) const
Computes the FFT or IFFT in-place using precomputed tables.
void apply_bluestein_transform(Array< Complex > &a, const bool invert, ThreadPool *pool, const size_t chunk_size) const
size_t size() const noexcept
Returns the transform size this plan was built for.
void apply_transform(Array< Complex > &a, const bool invert, ThreadPool *pool, const size_t chunk_size, const bool prefer_simd) const
bool should_use_neon(ThreadPool *pool, const bool prefer_simd) const noexcept
Array< Array< Real > > irfft_batch(const Array< Array< Complex > > &spectra) const
Computes a batch of Complex-to-Real Inverse FFTs.
Array< Array< Complex > > prfft_batch(ThreadPool &pool, const Array< Array< Real > > &input, const size_t chunk_size=0) const
Computes a parallel batch of Real-to-Complex FFTs.
Array< Complex > bluestein_kernel_forward_
void apply_mixed_radix_transform(Array< Complex > &a, const bool invert) const
std::shared_ptr< const Plan > bluestein_plan_
Array< Real > pinverse_transform_real(ThreadPool &pool, const Array< Complex > &input, const size_t chunk_size=0) const
Parallel IFFT projected back to real values.
Array< Complex > bluestein_kernel_inverse_
void apply_butterflies(Array< Complex > &a, const bool invert, ThreadPool *pool, const size_t chunk_size, const bool prefer_simd) const
Array< Complex > transformed(const Array< Complex > &input, const bool invert=false) const
Computes the FFT or IFFT returning a new array.
Array< Complex > inverse_transform(const Array< Complex > &input) const
Returns the inverse transform as a new array.
bool automatic_batch_simd_candidate(const size_t batch_size) const noexcept
Array< Array< Real > > pinverse_transform_real_batch(ThreadPool &pool, const Array< Array< Complex > > &input, const size_t chunk_size=0, const bool prefer_simd=true) const
Returns the parallel inverse batch transform projected back to real values.
bool supports_power_of_two_real_optimization() const noexcept
Array< Complex > pinverse_transform(ThreadPool &pool, const Array< Complex > &input, const size_t chunk_size=0) const
Parallel inverse transform returning a new array.
Array< Array< Real > > pirfft_batch(ThreadPool &pool, const Array< Array< Complex > > &spectra, const size_t chunk_size=0) const
Computes a parallel batch of Complex-to-Real Inverse FFTs.
Array< Complex > mixed_radix_transform_recursive(const Array< Complex > &input, const size_t factor_index, const size_t current_size, const bool invert) const
Array< Array< Complex > > inverse_transform_batch(const Array< Array< Complex > > &input) const
Returns the inverse batch transform as a new array of arrays.
bool should_use_avx2(ThreadPool *pool, const bool prefer_simd) const noexcept
Array< Complex > twiddles_
void ptransform_batch(ThreadPool &pool, Array< Array< Complex > > &batch, const bool invert, const size_t chunk_size=0, const bool prefer_simd=true) const
Parallel batch transform for equal-length complex inputs.
Array< Real > irfft(const Array< Complex > &spectrum) const
Computes the Complex-to-Real Inverse FFT (IRFFT).
Array< Array< Complex > > transformed_batch(const Array< Array< Complex > > &input, const bool invert=false, const bool prefer_simd=true) const
Returns the batch FFT/IFFT as a new array of arrays.
Array< Array< Complex > > pinverse_transform_batch(ThreadPool &pool, const Array< Array< Complex > > &input, const size_t chunk_size=0, const bool prefer_simd=true) const
Returns the parallel inverse batch transform.
Array< Array< Complex > > rfft_batch(const Array< Array< Real > > &input) const
Computes a batch of Real-to-Complex FFTs.
Array< Array< Real > > inverse_transform_real_batch(const Array< Array< Complex > > &input) const
Returns the inverse batch transform projected back to real values.
Array< Complex > bluestein_chirp_
void ptransform(ThreadPool &pool, Array< Complex > &a, const bool invert, const size_t chunk_size=0) const
Parallel in-place transform using a thread pool.
Plan()=default
Constructs an empty plan (size 0).
Stateful bank of SOS cascades, one per channel.
SOSFilter & channel(const size_t index)
SOSFilterBank(const size_t num_channels, const SectionsContainer §ions)
Array< Array< Real > > filter(const Array< Array< Real > > &signals)
Array< Real > filter_channel(const size_t index, const Container &signal)
const SOSFilter & channel(const size_t index) const
Array< Array< Real > > pfilter(ThreadPool &pool, const Array< Array< Real > > &signals, const size_t chunk_size=0)
Array< Real > filter_channel(const size_t index, const Array< Real > &signal)
Array< SOSFilter > filters_
size_t num_channels() const noexcept
void validate_channel_count(const size_t count, const char *ctx) const
Rationale: Internal safety guard.
SOSFilterBank(const size_t num_channels, const Array< BiquadSection > §ions)
Stateful cascade of second-order sections.
Array< Array< Real > > states_
size_t num_sections() const noexcept
Returns the number of sections in the cascade.
SOSFilter(const SectionsContainer §ions)
Construct a cascade of SOS from any biquad container.
SOSFilter(const Array< BiquadSection > §ions)
Construct a cascade of second-order sections (SOS).
bool configured() const noexcept
Returns whether the filter has been configured.
SOSFilter()=default
Default constructor for an unconfigured SOSFilter.
const Array< Real > & state(const size_t index) const
Returns the current delay-line state of section index.
Array< BiquadSection > sections_
Array< IIRCoefficients > coeffs_
Array< Real > filter(const Array< Real > &signal)
Processes a signal block through the stateful SOS cascade.
Array< Real > filter(const Container &signal)
Processes a real-valued container through the stateful SOS cascade.
void reset()
Resets all internal delay lines to zero.
void require_configured(const char *ctx) const
Stateful STFT processor for chunked real-time analysis.
Array< Array< Complex > > pflush(ThreadPool &pool, const size_t chunk_size=0)
size_t frame_size() const noexcept
void initialize_pending()
void append_samples(const Array< Real > &block)
size_t fft_size() const noexcept
bool centered() const noexcept
STFTProcessor(const Array< Real > &window, const STFTOptions &options)
Construct an STFT processor with a fixed analysis window.
bool finalized() const noexcept
Array< Array< Complex > > pprocess_block(ThreadPool &pool, const Container &block, const size_t chunk_size=0)
Array< Array< Complex > > flush()
Array< Array< Complex > > process_block(const Container &block)
void require_configured(const char *ctx) const
STFTProcessor(const size_t frame_size, const STFTOptions &options)
Construct an STFT processor using a Hann window.
Array< Array< Complex > > process_block(const Array< Real > &block)
const Array< Real > & window() const noexcept
Array< Array< Complex > > emit_ready_frames(ThreadPool *pool, const size_t chunk_size, const bool allow_partial_frames)
Array< Array< Complex > > pprocess_block(ThreadPool &pool, const Array< Real > &block, const size_t chunk_size=0)
STFTProcessor(const WindowContainer &window, const STFTOptions &options)
Construct an STFT processor from any real container.
size_t hop_size() const noexcept
STFTProcessor()=default
Default constructor for an unconfigured processor.
bool pad_end() const noexcept
bool configured() const noexcept
Fast Fourier Transform (FFT) and DSP Toolkit.
static Array< Array< Array< Complex > > > inverse_transform3d(const Array< Array< Array< Complex > > > &input)
Functional inverse 3-D FFT wrapper.
static Array< Array< Complex > > rfft_batch(const Array< Array< Real > > &input)
Functional compact real FFT for equal-length real batches.
static Array< Real > group_delay(const SectionsContainer §ions, const size_t num_points=512, const bool whole=false)
static Array< Complex > spectrum(const Container &input)
DSP alias for complex-valued container forward FFT.
static Array< Real > istft(const Array< Array< Complex > > &spectrogram, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options)
Reconstructs a real signal using explicit ISTFT options.
static Array< Real > signed_binomial(const size_t power, const Real sign)
Rationale: Computes (1 + sign*x)^power coefficients.
static Array< Complex > transform_padded(const Container &input)
Forward FFT with padding for real-valued containers.
static void validate_real_spectrum(const Array< Complex > &input, const char *ctx)
Rationale: Ensures a complex spectrum has the Hermitian symmetry (X[k] = X[N-k]*) required to be the ...
static Array< BiquadSection > bessel_bandstop(const size_t order, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Bessel band-stop design returned as SOS.
static bool has_near_pole_zero_cancellation(const Array< Real > &numerator, const Array< Real > &denominator, const Real tolerance)
static Array< Array< Complex > > stft(const Array< Real > &signal, const Array< Real > &window, const size_t hop_size, const bool pad_end=true)
Computes a basic STFT for a real signal using a custom window.
static WeightedFrequencyGrid build_weighted_frequency_grid(const size_t num_taps, const Array< Real > &bands, const Array< Real > &desired, const Real sample_rate, const Array< Real > &weights, const size_t grid_density, const char *ctx)
static Real hermitian_tolerance(const Complex &lhs, const Complex &rhs, const size_t n) noexcept
Rationale: Numerical tolerance for Hermitian symmetry checks.
static Array< Array< Real > > frame_signal(const Array< Real > &signal, const size_t frame_size, const size_t hop_size, const bool pad_end=true)
Splits a real signal into frames using a fixed hop size.
static Array< Real > apply_blackman_window(const Array< Real > &signal)
Applies a Blackman window of matching size to a real signal.
static Array< Real > phase_delay(const Array< Real > &numerator, const size_t num_points=512, const bool whole=false)
static Array< Array< Real > > pbatched_sosfilt(ThreadPool &pool, const Array< Array< Real > > &signals, const Array< BiquadSection > §ions, const size_t chunk_size=0)
and static Is_Real_Container< Container2 > Array< Real > overlap_save_convolution(const Container1 &signal, const Container2 &kernel, const size_t block_size=0)
Overlap-save convolution for real-valued containers.
static Array< Real > trim_leading_zeros_copy(const Array< Real > &input) noexcept
Rationale: Copy utility that removes leading zero coefficients, effectively normalizing the polynomia...
static Array< Real > hamming_window(const size_t n)
Returns a Hamming window of length n.
and Is_Real_Container< DesiredContainer > and static Is_Real_Container< WeightContainer > Array< Real > remez(const size_t num_taps, const BandContainer &bands, const DesiredContainer &desired, const Real sample_rate, const WeightContainer &weights, const size_t grid_density=32, const size_t max_iterations=64)
SimdBackend
SIMD hardware acceleration backends.
@ scalar
Generic implementation.
static Array< size_t > frame_offsets(const size_t signal_size, const size_t frame_size, const size_t hop_size, const bool pad_end=true)
Returns the frame start offsets used by frame_signal.
static std::pair< Real, Real > prewarp_band_edges(const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate, const char *ctx)
Prewarps a pair of band edges for bilinear transform.
static Array< Real > sos_filtfilt_impl(const Array< Real > &signal, const Array< BiquadSection > §ions, const char *ctx)
static void validate_no_near_pole_zero_cancellation(const Array< BiquadSection > §ions, const Real tolerance)
static Array< Array< Complex > > compact_real_batch_spectra(const Array< Array< Complex > > &full_spectra)
static Array< Real > project_real_output(const Array< Complex > &input, const char *ctx)
Projects a complex-valued output array back to the real domain.
static Array< Complex > build_complex_input(const Container &input)
static Array< Real > phase_delay(const BiquadSection §ion, const size_t num_points=512, const bool whole=false)
static Array< Complex > ptransformed(ThreadPool &pool, const Array< Complex > &input, const bool invert=false, const size_t chunk_size=0)
Parallel FFT/IFFT that returns a new array.
and static Is_Real_Container< DenContainer > GainMarginInfo gain_margin(const NumContainer &numerator, const DenContainer &denominator, const size_t num_points=1024, const bool whole=false)
static Array< Real > lfilter(const Array< Real > &signal, const Array< Real > &numerator, const Array< Real > &denominator, const Array< Real > &initial_state={})
One-shot causal IIR filtering.
and static Is_Real_Container< DenContainer > Real minimum_pole_zero_distance(const NumContainer &numerator, const DenContainer &denominator)
static constexpr const char * simd_preference_name(const SimdPreference preference) noexcept
Returns the human-readable name of a SIMD preference.
static Array< Array< Complex > > stft(const Array< Real > &signal, const Array< Real > &window, const STFTOptions &options)
Computes an STFT with explicit analysis options.
static FrequencyResponse freqz_impl(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points, const bool whole, const char *ctx)
static Array< Real > apply_blackman_window(const Container &signal)
SimdPreference
User-selected SIMD dispatch preference.
@ scalar_only
Force the portable scalar implementation.
@ automatic
Use the fastest available hardware.
@ neon_only
Force NEON (must be supported by CPU).
@ avx2_only
Force AVX2 (must be supported by CPU).
static Array< Real > phase_delay_impl(const FrequencyResponse &response)
implementation of numeric phase delay calculation.
static PhaseMarginInfo phase_margin(const BiquadSection §ion, const size_t num_points=1024, const bool whole=false)
and static Is_Real_Container< Container2 > Array< Real > multiply(const Container1 &a, const Container2 &b)
Multiplication for real-valued containers.
static Real modified_bessel_i0(const Real x) noexcept
Standard Zeroth-order Modified Bessel function of the first kind.
static Array< Complex > rfft(const Container &input)
Compact real FFT for real-valued containers.
static constexpr size_t saturating_product(const size_t lhs, const size_t rhs) noexcept
static Array< BiquadSection > butterworth_highpass(const size_t order, const Real cutoff_frequency, const Real sample_rate)
Digital Butterworth high-pass design returned as SOS.
static PhaseMarginInfo phase_margin(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points=1024, const bool whole=false)
static size_t resolve_irfft_signal_size(const Array< Complex > &spectrum, const size_t signal_size, const char *ctx)
static Array< Real > filtfilt(const Array< Real > &signal, const BiquadSection §ion)
Zero-phase filtering of one biquad section.
static Array< Real > analytic_phase_delay_impl(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points, const bool whole, const char *ctx)
implementation of analytic phase delay calculation.
static Real evaluate_cosine_series(const Array< Real > &coefficients, const Real omega) noexcept
static Array< BiquadSection > bessel_lowpass(const size_t order, const Real cutoff_frequency, const Real sample_rate)
Digital Bessel low-pass design returned as SOS.
static Array< Complex > pspectrum(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of spectrum(const Container&).
static Array< Complex > poles(const IIRCoefficients &coeffs)
static void validate_stable(const Array< Real > &denominator, const Real min_margin)
and static Is_Real_Container< WindowContainer > PowerSpectralDensity welch(const SignalContainer &signal, const WindowContainer &window, const Real sample_rate, const WelchOptions &options={})
static Array< Complex > apply_hann_window(const Array< Complex > &signal)
Applies a Hann window of matching size to a complex signal.
and Is_Real_Container< ContainerY > and static Is_Real_Container< WindowContainer > CoherenceEstimate coherence(const ContainerX &x, const ContainerY &y, const WindowContainer &window, const Real sample_rate, const WelchOptions &options={})
static Array< Real > solve_dense_system(Array< Real > matrix, Array< Real > rhs, const size_t n, const char *ctx)
static Array< Real > filtfilt(const Array< Real > &signal, const IIRCoefficients &coeffs)
Zero-phase IIR filtering with explicit coefficient groups.
static void validate_no_near_pole_zero_cancellation(const IIRCoefficients &coeffs, const Real tolerance)
static BiquadSection section_from_coefficients(const Array< Real > &numerator, const Array< Real > &denominator, const char *ctx)
Rationale: Converts generic IIR coefficients into a stabilized BiquadSection structure,...
static bool roots_pass_residual_check(const Array< Real > &coefficients, const Array< Complex > &roots, const Real tol) noexcept
static Array< Real > firwin_bandstop(const size_t num_taps, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate, const Array< Real > &window)
FIR band-stop design via the window method.
static Array< Real > solve_remez_cosine_series(const Array< Real > &omega, const Array< Real > &desired, const Array< Real > &weight, const Array< size_t > &extrema, const size_t half_order, const char *ctx)
static Array< Real > pistft(ThreadPool &pool, const Array< Array< Complex > > &spectrogram, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size, const size_t signal_length=0, const size_t chunk_size=0)
Parallel STFT inversion using a thread pool.
static Array< Array< Real > > overlap_add_convolution_batch(const Array< Array< Real > > &signals, const Array< Real > &kernel, const size_t block_size=0)
Multichannel overlap-add convolution with a shared kernel FFT.
static void validate_no_near_pole_zero_cancellation(const BiquadSection §ion, const Real tolerance)
static Array< Complex > apply_window(const Array< Complex > &signal, const Array< Real > &window)
Applies a real window sample-by-sample to a complex signal.
static Real stability_margin(const SectionsContainer §ions)
static Array< Real > lfilter(const Array< Real > &signal, const IIRCoefficients &coeffs, const Array< Real > &initial_state={})
static bool has_near_pole_zero_cancellation(const Array< Complex > &zeros, const Array< Complex > &poles, const Real tolerance)
Detects near pole/zero cancellations under a tolerance.
static Array< Real > overlap_save_convolution(const Array< Real > &signal, const Array< Real > &kernel, const size_t block_size=0)
Convenience wrapper for overlap-save real convolution.
static Array< Complex > zeros(const BiquadSection §ion)
static Array< Real > phase_spectrum(const Container &input)
static constexpr bool is_power_of_two(const size_t n) noexcept
Checks if a given number is a power of two.
static Array< BiquadSection > chebyshev2_lowpass(const size_t order, const Real attenuation_db, const Real cutoff_frequency, const Real sample_rate)
Digital Chebyshev-II low-pass design returned as SOS.
static void validate_stable(const SectionsContainer §ions)
static void append_all(Array< T > &dst, const Array< T > &src)
static SeriesEvaluation evaluate_series_at_unit_circle(const Array< Real > &coefficients, const Real omega, const char *ctx)
static void ptransform(ThreadPool &pool, Array< Complex > &a, const bool invert, const size_t chunk_size=0)
Parallel in-place FFT/IFFT using a thread pool.
static Array< Array< Complex > > pstft(ThreadPool &pool, const Array< Real > &signal, const size_t frame_size, const STFTOptions &options, const size_t chunk_size=0)
Parallel Hann-window STFT with explicit analysis options.
static Array< Real > firwin_bandpass(const size_t num_taps, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate, const Real attenuation_db)
FIR band-pass design using a Kaiser window.
static bool try_durand_kerner_roots(const Array< Real > &coefficients, const Real tol, Array< Complex > &roots, const size_t max_iterations=192) noexcept
static bool try_aberth_ehrlich_roots(const Array< Real > &coefficients, const Real tol, Array< Complex > &roots, const size_t max_iterations=192) noexcept
static Array< Complex > zeros(const Array< Real > &numerator)
Returns the finite zeros of a transfer numerator.
static Array< Array< Array< Complex > > > batched_stft(const Array< Array< Real > > &signals, const size_t frame_size, const STFTOptions &options)
Hann-window batched STFT.
and static Is_Real_Container< WindowContainer > Array< Array< Complex > > stft(const SignalContainer &signal, const WindowContainer &window, const STFTOptions &options)
static IIRCoefficients apply_analog_rational_transform(const Array< Real > &numerator, const Array< Real > &denominator, const Array< Real > &map_numerator, const Array< Real > &map_denominator, const char *ctx)
Rationale: Transforms an analog transfer function H(s) via rational substitution s = N(z)/D(z),...
static Real refine_scalar_crossing(const Evaluator &evaluator, const Real x0, const Real y0, const Real x1, const Real y1, const Real target)
static Array< BiquadSection > transfer_function_to_sections(const IIRCoefficients &coeffs, const Array< RootGroup > &zero_groups, const char *ctx)
Rationale: Decomposes a general IIR transfer function into a cascade of Second-Order Sections (SOS) b...
static Array< Real > multiply_real_optimized(const Array< Real > &a, const Array< Real > &b, ThreadPool *pool=nullptr, const size_t chunk_size=0)
Rationale: Efficiently multiplies two real sequences using a single N-point complex FFT by packing on...
static GainMarginInfo gain_margin(const FrequencyResponse &response)
Estimates gain margin around the -pi phase crossover.
static Real comp_ellint_1_impl(const Real k)
static Array< Array< Complex > > stft(const Array< Real > &signal, const size_t frame_size, const size_t hop_size, const bool pad_end=true)
Computes a basic STFT for a real signal using a Hann window.
static Array< Real > phase_delay(const FrequencyResponse &response)
Estimates phase delay from a sampled frequency response.
static void validate_overlap_constraints(const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size, const bool validate_nola, const bool validate_cola, const char *ctx)
Rationale: Internal validator for STFT/ISTFT windowing parameters.
static Real stability_margin(const Array< Real > &denominator)
Signed stability margin relative to the unit circle.
static Array< Real > firwin_lowpass(const size_t num_taps, const Real cutoff_frequency, const Real sample_rate, const Real attenuation_db)
FIR low-pass design via the window method using a Kaiser window.
static Array< Real > istft_impl(const Array< Array< Complex > > &spectrogram, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options, ThreadPool *pool, const size_t chunk_size)
implementation of the Inverse Short-Time Fourier Transform.
static Real minimum_pole_zero_distance(const Array< Complex > &zeros, const Array< Complex > &poles)
Returns the minimum pole/zero distance in a greedy pairing.
static Array< Complex > zeros(const IIRCoefficients &coeffs)
static Real window_enbw(const Container &window)
static Real align_phase_near_reference(const Real raw_phase, const Real reference_phase) noexcept
static Array< Real > pistft(ThreadPool &pool, const Array< Array< Complex > > &spectrogram, const size_t frame_size, const size_t hop_size, const size_t signal_length=0, const size_t chunk_size=0)
Parallel STFT inversion using a Hann window pair.
static Array< Complex > apply_hamming_window(const Array< Complex > &signal)
Applies a Hamming window of matching size to a complex signal.
static PhaseMarginInfo phase_margin_impl(const FrequencyResponse &response)
static Array< Real > inverse_transform_real_general_impl(const Array< Complex > &input, const char *ctx, ThreadPool *pool, const size_t chunk_size)
static Array< Real > build_real_input(const Container &input)
static Array< Real > firls_impl(const size_t num_taps, const Array< Real > &bands, const Array< Real > &desired, const Real sample_rate, const Array< Real > &weights, const char *ctx)
implementation of Least-Squares FIR design.
static Array< Real > filtfilt(const Array< Real > &signal, const Array< BiquadSection > §ions)
Zero-phase filtering of a cascade of biquad sections.
and static Is_Real_Container< Container2 > Array< Real > poverlap_add_convolution(ThreadPool &pool, const Container1 &signal, const Container2 &kernel, const size_t block_size=0, const size_t chunk_size=0)
Parallel version of overlap_add_convolution(const Container1&, const Container2&, size_t).
static FrequencyResponse freqz(const BiquadSection §ion, const size_t num_points=512, const bool whole=false)
static Array< Real > analytic_group_delay_impl(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points, const bool whole, const char *ctx)
Rationale: Analytic calculation of group delay (-dphi/domega) using the complex derivative of the tra...
static Array< Real > pirfft(ThreadPool &pool, const Array< Complex > &spectrum, const size_t signal_size=0, const size_t chunk_size=0)
Parallel version of irfft(const Array<Complex>&, size_t).
static Array< Complex > zeros(const Container &numerator)
static CoherenceEstimate coherence(const Array< Real > &x, const Array< Real > &y, const Array< Real > &window, const Real sample_rate, const WelchOptions &options={})
Magnitude-squared coherence estimate using Welch averages.
static Array< Real > unwrap_phase(const Array< Real > &phase)
implementation of phase unwrapping (removes +/- 2*pi jumps).
static Array< Real > polynomial_power(const Array< Real > &poly, const size_t exponent, const char *ctx)
implementation of P(x)^n via successive multiplication.
static bool is_stable(const Container &denominator)
static constexpr bool Is_Real_Batch_Container
static Array< Real > phase_delay(const Array< BiquadSection > §ions, const size_t num_points=512, const bool whole=false)
static bool neon_dispatch_available() noexcept
Returns whether the runtime CPU can execute the NEON kernel when it has been compiled in.
static Array< Real > build_padded_frame_from_prefix(const Array< Real > &input, const size_t frame_size)
static Array< Array< Complex > > ptransformed2d(ThreadPool &pool, const Array< Array< Complex > > &input, const bool invert=false, const size_t chunk_size=0)
Parallel functional 2-D FFT/IFFT wrapper for matrices.
static void scatter_axis_slice(Array< Complex > &data, const size_t base_offset, const size_t axis_stride, const Array< Complex > &slice)
Rationale: Writes back a transformed contiguous slice into its original (potentially non-contiguous) ...
static Array< Array< Real > > pbatched_lfilter(ThreadPool &pool, const Array< Array< Real > > &signals, const Array< Real > &numerator, const Array< Real > &denominator, const Array< Array< Real > > &initial_states={}, const size_t chunk_size=0)
static constexpr const char * simd_backend_name(const SimdBackend backend) noexcept
Returns the human-readable name of a SIMD backend.
and static Is_Real_Container< WindowContainer > Array< Array< Complex > > stft(const SignalContainer &signal, const WindowContainer &window, const size_t hop_size, const bool pad_end=true)
static Array< Complex > poles(const Array< BiquadSection > §ions)
static Array< Array< Real > > pbatched_istft(ThreadPool &pool, const Array< Array< Array< Complex > > > &spectrograms, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options, const Array< size_t > &signal_lengths={}, const size_t chunk_size=0)
Parallel batched ISTFT across signals.
static Array< Array< Complex > > pinverse_transform2d(ThreadPool &pool, const Array< Array< Complex > > &input, const size_t chunk_size=0)
Parallel functional inverse 2-D FFT wrapper.
static Array< BiquadSection > chebyshev2_highpass(const size_t order, const Real attenuation_db, const Real cutoff_frequency, const Real sample_rate)
Digital Chebyshev-II high-pass design returned as SOS.
static Array< Real > irfft(const Array< Complex > &spectrum, const size_t signal_size=0)
Reconstructs a real signal from a compact rfft() spectrum.
static Array< BiquadSection > design_transformed_sections(const AnalogPrototype &prototype, const Array< Real > &map_numerator, const Array< Real > &map_denominator, const Real sample_rate, const char *ctx)
Generic implementation for designing transformed SOS cascades.
static size_t resolve_welch_fft_size(const WelchOptions &options, const size_t frame_size, const char *ctx)
Default FFT size for Welch analysis (zero-padded to power of 2).
static Array< BiquadSection > butterworth_lowpass(const size_t order, const Real cutoff_frequency, const Real sample_rate)
Digital Butterworth low-pass design returned as SOS.
static Array< Array< Array< Complex > > > pbatched_stft(ThreadPool &pool, const Array< Array< Real > > &signals, const size_t frame_size, const STFTOptions &options, const size_t chunk_size=0)
Parallel Hann-window batched STFT.
static bool is_stable(const Array< Real > &denominator)
Checks BIBO stability from denominator roots.
static Array< Array< Complex > > prfft_batch(ThreadPool &pool, const Array< Array< Real > > &input, const size_t chunk_size=0)
Parallel compact real FFT for equal-length real batches.
static Array< Real > group_delay(const Array< Real > &numerator, const size_t num_points=512, const bool whole=false)
static BalancedPolynomial balance_polynomial_for_roots(const Array< Real > &coefficients, const char *ctx)
static Array< Complex > polynomial_roots_impl(const Array< Real > &coefficients, const char *ctx)
Rationale: Master implementation for polynomial root-finding.
static Real minimum_cancellation_distance_impl(const Array< PoleZeroPair > &pairs) noexcept
Rationale: Estimates the minimum distance between any matched pole and zero to identify potential can...
static Array< Real > phase_delay(const SectionsContainer §ions, const size_t num_points=512, const bool whole=false)
static Array< Array< Real > > pmultichannel_istft(ThreadPool &pool, const Array< Array< Array< Complex > > > &spectrograms, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options={}, const Array< size_t > &signal_lengths={}, const SpectrogramLayout layout=SpectrogramLayout::channel_frame_bin, const size_t chunk_size=0)
Parallel multichannel ISTFT from either supported layout.
static TensorLayout row_major_layout(const Array< size_t > &shape)
Builds a default row-major tensor layout for a flat buffer.
static constexpr size_t twiddle_refresh_period
static Array< Real > hann_window(const size_t n)
Returns a Hann window of length n.
static Array< Complex > rfft(const Array< Real > &input)
Computes the compact real FFT of size floor(N/2) + 1.
static Real integrate_offset_cos_basis(const Real omega_lo, const Real omega_hi, const size_t harmonic) noexcept
implementation helper for firls slope integration.
static Array< Real > remez(const size_t num_taps, const Array< Real > &bands, const Array< Real > &desired, const Real sample_rate, const Array< Real > &weights={}, const size_t grid_density=32, const size_t max_iterations=64)
FIR equiripple design via a dense-grid Remez exchange.
static Array< Array< Array< Complex > > > pinverse_transform3d(ThreadPool &pool, const Array< Array< Array< Complex > > > &input, const size_t chunk_size=0)
Parallel functional inverse 3-D FFT wrapper.
static Array< Complex > pspectrum(ThreadPool &pool, const Array< Complex > &input, const size_t chunk_size=0)
Parallel version of spectrum(const Array<Complex>&).
static Array< Real > pistft(ThreadPool &pool, const Array< Array< Complex > > &spectrogram, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options, const size_t chunk_size=0)
Parallel ISTFT with explicit options.
static Real minimum_pole_zero_distance(const Array< BiquadSection > §ions)
static Array< Real > analytic_sos_group_delay_impl(const Array< BiquadSection > §ions, const FrequencyResponse &response, const char *ctx)
implementation of analytic group delay for SOS cascades.
static Array< size_t > evenly_spaced_extrema(const size_t grid_size, const size_t extremal_count)
static Array< Real > trim_polynomial_leading_zeros(const Array< Real > &input, const char *ctx)
Validation wrapper for trim_leading_zeros_copy.
static Array< Real > power_spectrum(const Container &input)
static Real stability_margin(const Array< BiquadSection > §ions)
static Array< Complex > pinverse_transform(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of inverse_transform(const Container&).
and static Is_Real_Container< CoeffContainer > Array< Real > resample_poly(const SignalContainer &signal, const size_t up, const size_t down, const CoeffContainer &coeffs)
static Array< Array< Real > > project_real_batch_output(const Array< Array< Complex > > &input, const char *ctx)
static Real evaluate_fir_response_magnitude(const Array< Real > &coeffs, const Real omega, const char *ctx)
implementation helper for FIR frequency response magnitude.
and static Is_Real_Container< DesiredContainer > Array< Real > remez(const size_t num_taps, const BandContainer &bands, const DesiredContainer &desired, const Real sample_rate)
static Array< Real > irfft(const Container &spectrum, const size_t signal_size=0)
Inverse real transform (compact) for complex-valued containers.
static Array< Array< Array< Complex > > > ptransformed2d_batch(ThreadPool &pool, const Array< Array< Array< Complex > > > &input, const bool invert=false, const size_t chunk_size=0)
Parallel functional 2-D batched FFT/IFFT wrapper.
and static Is_Real_Container< Container2 > Array< Real > pmultiply(ThreadPool &pool, const Container1 &a, const Container2 &b, const size_t chunk_size=0)
Parallel version of multiply(const Container1&, const Container2&).
static AnalogPrototype butterworth_prototype(const size_t order, const char *ctx)
static Array< Real > overlap_profile(const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size, const char *ctx)
Rationale: Computes the effective window product for a windowed overlap-add process to validate recon...
static Array< Array< Real > > batched_istft(const Array< Array< Array< Complex > > > &spectrograms, const Array< Real > &window, const ISTFTOptions &options, const Array< size_t > &signal_lengths={})
Batched ISTFT using a shared analysis/synthesis window.
static Real max_root_radius(const Array< Complex > &roots) noexcept
Returns the magnitude of the largest root in the set.
static Array< Array< Array< Complex > > > ptransformed3d(ThreadPool &pool, const Array< Array< Array< Complex > > > &input, const bool invert=false, const size_t chunk_size=0)
Parallel functional 3-D FFT/IFFT wrapper for tensors.
static Array< Complex > prfft(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of rfft(const Container&).
static SimdBackend simd_backend() noexcept
Returns the default SIMD backend used by standalone plan transforms for this precision.
static Array< Array< Real > > poverlap_add_convolution_batch(ThreadPool &pool, const Array< Array< Real > > &signals, const Array< Real > &kernel, const size_t block_size=0, const size_t chunk_size=0)
Parallel version of overlap_add_convolution_batch.
static Real normalized_sinc(const Real x) noexcept
Standard normalized sinc function sin(pi*x)/(pi*x).
static Real ellint_1_impl(const Real k, const Real phi)
static Array< Real > demean_copy(const Array< Real > &input)
Returns a copy with DC offset removed.
static Array< Array< Real > > frame_signal(const Container &signal, const size_t frame_size, const size_t hop_size, const bool pad_end=true)
static const char * detected_simd_backend_name() noexcept
Returns the detected hardware SIMD backend name.
static Array< Array< Complex > > expand_real_batch_spectra(const Array< Array< Complex > > &spectra, const size_t signal_size, const char *ctx)
static Array< Real > inverse_transform_real_optimized_impl(const Array< Complex > &input, const char *ctx, ThreadPool *pool, const size_t chunk_size, InversePackedTransform inverse_packed)
Rationale: Optimized Inverse FFT for real signals of size N (power-of-two).
static Array< Real > resample_poly(const SignalContainer &signal, const size_t up, const size_t down, const ResamplePolyOptions &options={})
static Array< Real > lfilter(const Array< Real > &signal, const BiquadSection §ion, const Array< Real > &initial_state={})
static Array< Real > istft(const Array< Array< Complex > > &spectrogram, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size, const size_t signal_length=0)
Reconstructs a real signal from an STFT using overlap-add.
static Array< PoleZeroPair > pair_poles_and_zeros(const Array< BiquadSection > §ions)
static Array< Real > pinverse_transform_real(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of inverse_transform_real(const Container&).
static Real stability_margin(const Container &denominator)
static CrossSpectralDensity csd(const Array< Real > &x, const Array< Real > &y, const Array< Real > &window, const Real sample_rate, const WelchOptions &options={})
One-sided cross-spectral density estimate using a custom window.
static Array< Complex > zeros(const SectionsContainer §ions)
static Array< Array< Real > > prepare_welch_frames(const Array< Real > &signal, const Array< Real > &window, const WelchOptions &options, const char *ctx)
Rationale: Encapsulates frame extraction, detrending and windowing common to Welch and other spectral...
static Complex evaluate_transfer_at(const Array< Real > &numerator, const Array< Real > &denominator, const Real omega, const char *ctx)
static Array< BiquadSection > elliptic_bandstop(const size_t order, const Real ripple_db, const Real attenuation_db, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital elliptic/Cauer band-stop design returned as SOS.
static bool is_stable(const SectionsContainer §ions)
static Array< Real > one_sided_frequency_grid(const size_t fft_size, const Real sample_rate, const char *ctx)
Maps bin indices to frequency in Hz.
static Array< Complex > lift_real_input(const Array< Real > &input)
Lifts a real-valued input array to the complex domain.
static Array< Real > pfiltfilt(ThreadPool &pool, const Array< Real > &signal, const Array< Real > &coeffs, const size_t block_size=0, const size_t chunk_size=0)
Parallel zero-phase FIR filtering via forward-backward convolution.
static constexpr bool Is_Biquad_Container
static Array< Array< Real > > batched_istft(const Array< Array< Array< Complex > > > &spectrograms, const size_t frame_size, const ISTFTOptions &options, const Array< size_t > &signal_lengths={})
Hann-window batched ISTFT.
static Array< PoleZeroPair > pair_poles_and_zeros(const Array< Real > &numerator, const Array< Real > &denominator)
static Array< Array< Real > > pbatched_istft(ThreadPool &pool, const Array< Array< Array< Complex > > > &spectrograms, const size_t frame_size, const ISTFTOptions &options, const Array< size_t > &signal_lengths={}, const size_t chunk_size=0)
Parallel Hann-window batched ISTFT.
static size_t validate_stft_options(const Array< Real > &window, const STFTOptions &options, const char *ctx)
static CoherenceEstimate coherence(const Array< Real > &x, const Array< Real > &y, const size_t frame_size, const Real sample_rate, const WelchOptions &options={})
Magnitude-squared coherence estimate using a Hann window.
static Array< BiquadSection > elliptic_highpass(const size_t order, const Real ripple_db, const Real attenuation_db, const Real cutoff_frequency, const Real sample_rate)
Digital elliptic/Cauer high-pass design returned as SOS.
static void bit_reverse(Array< Complex > &a) noexcept
Performs bit-reversal permutation in-place.
static Complex evaluate_sos_transfer_at(const Array< BiquadSection > §ions, const Real omega, const char *ctx)
static Array< BiquadSection > bessel_highpass(const size_t order, const Real cutoff_frequency, const Real sample_rate)
Digital Bessel high-pass design returned as SOS.
static Array< Array< Complex > > pstft(ThreadPool &pool, const Container &signal, const size_t frame_size, const STFTOptions &options, const size_t chunk_size=0)
static GainMarginInfo gain_margin(const SectionsContainer §ions, const size_t num_points=1024, const bool whole=false)
static Array< Real > apply_hann_window(const Array< Real > &signal)
Applies a Hann window of matching size to a real signal.
static Array< Array< Array< Complex > > > transpose_spectrogram_layout_impl(const Array< Array< Array< Complex > > > &input, const SpectrogramLayout source, const SpectrogramLayout target, const char *ctx)
Rationale: Internal implementation for transposing between multichannel spectrogram layouts (e....
static Array< Array< Complex > > pstft(ThreadPool &pool, const Array< Real > &signal, const Array< Real > &window, const STFTOptions &options, const size_t chunk_size=0)
Parallel STFT using explicit analysis options.
static void normalize_fir_at_omega(Array< Real > &coeffs, const Real omega, const char *ctx)
Normalizes FIR coefficients to have unit gain at frequency omega.
static Array< Complex > spectrum(const Array< Complex > &input)
DSP alias for the forward FFT.
static Array< Real > group_delay(const Array< BiquadSection > §ions, const size_t num_points=512, const bool whole=false)
static Real stability_margin(const IIRCoefficients &coeffs)
static Real interpolate_crossing(const Real x0, const Real y0, const Real x1, const Real y1, const Real target) noexcept
static Array< Real > pirfft(ThreadPool &pool, const Container &spectrum, const size_t signal_size=0, const size_t chunk_size=0)
Parallel version of irfft(const Container&, size_t).
static Array< Real > apply_hamming_window(const Container &signal)
static Array< Array< Complex > > stft(const Container &signal, const size_t frame_size, const STFTOptions &options)
static Array< Real > sosfilt(const Array< Real > &signal, const Array< BiquadSection > §ions)
One-shot causal filtering of a cascade of second-order sections.
static Array< Complex > transform(const Container &input)
Forward FFT for complex-valued containers.
static Array< Real > upfirdn(const Array< Real > &signal, const Array< Real > &coeffs, const size_t up=1, const size_t down=1)
Polyphase upsample-filter-downsample for real signals.
static void transform_axes_impl(Array< Complex > &data, const TensorLayout &layout, const Array< size_t > &axes, const bool invert, ThreadPool *pool=nullptr, const size_t chunk_size=0)
Rationale: Internal implementation of a multi-axis tensor transform.
static SimdBackend detected_simd_backend() noexcept
Returns the best SIMD backend supported by the current hardware among the kernels compiled for this p...
static Array< BiquadSection > butterworth_bandpass(const size_t order, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Butterworth band-pass design returned as SOS.
static Array< Complex > inverse_transform(const Container &input)
Forward IFFT for complex-valued containers.
static Array< Complex > ptransformed_axes(ThreadPool &pool, const Array< Complex > &input, const TensorLayout &layout, const Array< size_t > &axes, const bool invert=false, const size_t chunk_size=0)
Parallel functional tensor FFT/IFFT wrapper for flat buffers.
static Array< T > prefix_copy(const Array< T > &input, const size_t length)
Functional prefix extraction.
static Array< Array< Real > > batched_lfilter(const Array< Array< Real > > &signals, const Array< Real > &numerator, const Array< Real > &denominator, const Array< Array< Real > > &initial_states={})
One-shot batched causal IIR filtering across real channels.
and static Is_Real_Container< WindowContainer > Array< Complex > windowed_spectrum(const SignalContainer &signal, const WindowContainer &window)
static Real window_energy(const Container &window)
static Array< Array< Real > > multichannel_istft(const Array< Array< Array< Complex > > > &spectrograms, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options={}, const Array< size_t > &signal_lengths={}, const SpectrogramLayout layout=SpectrogramLayout::channel_frame_bin)
Multichannel ISTFT from either channel-major or frame-major spectrogram layouts.
static constexpr bool avx2_kernel_compiled() noexcept
Returns whether the AVX2 kernel was compiled for this precision.
static Array< Real > group_delay(const FrequencyResponse &response)
Estimates group delay from a sampled frequency response.
static FrequencyResponse freqz(const Array< BiquadSection > §ions, const size_t num_points=512, const bool whole=false)
static Real interpolate_value_at(const Real x0, const Real y0, const Real x1, const Real y1, const Real x) noexcept
static Complex group_center(const Array< Complex > &roots) noexcept
Geometric center of a root set.
static Array< Real > istft(const Array< Array< Complex > > &spectrogram, const Array< Real > &window, const ISTFTOptions &options)
Reconstructs a real signal using one shared analysis/synthesis window.
and Is_Real_Container< ContainerY > and static Is_Real_Container< WindowContainer > CrossSpectralDensity csd(const ContainerX &x, const ContainerY &y, const WindowContainer &window, const Real sample_rate, const WelchOptions &options={})
static Array< Complex > apply_blackman_window(const Array< Complex > &signal)
Applies a Blackman window of matching size to a complex signal.
static Array< Real > firwin_highpass(const size_t num_taps, const Real cutoff_frequency, const Real sample_rate, const Real attenuation_db)
FIR high-pass design using a Kaiser window.
static void validate_stable(const Array< BiquadSection > §ions)
static constexpr bool neon_kernel_compiled() noexcept
Returns whether the NEON kernel was compiled for this precision.
static FrequencyResponse freqz(const IIRCoefficients &coeffs, const size_t num_points=512, const bool whole=false)
static void validate_stable(const IIRCoefficients &coeffs, const Real min_margin)
static Array< BiquadSection > elliptic_bandpass(const size_t order, const Real ripple_db, const Real attenuation_db, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital elliptic/Cauer band-pass design returned as SOS.
static GainMarginInfo gain_margin(const BiquadSection §ion, const size_t num_points=1024, const bool whole=false)
and static Is_Real_Container< CoeffContainer > Array< Real > pfiltfilt(ThreadPool &pool, const SignalContainer &signal, const CoeffContainer &coeffs, const size_t block_size=0, const size_t chunk_size=0)
static bool avx2_runtime_available() noexcept
Returns whether AVX2 dispatch is supported at runtime.
static bool avx2_dispatch_available() noexcept
Returns whether the runtime CPU can execute the AVX2 double kernel when it has been compiled in.
static SimdBackend batched_plan_simd_backend() noexcept
Returns the SIMD backend selected for throughput-oriented batch plan paths under the current runtime ...
static Array< Real > filtfilt(const SignalContainer &signal, const BiquadSection §ion)
static Real window_coherent_gain(const Container &window)
static Array< Complex > multiply(const Array< Complex > &a, const Array< Complex > &b)
Multiplies two complex-valued sequences using convolution.
static std::pair< Array< Real >, Real > divide_polynomial_by_linear_root(const Array< Real > &coefficients, const Real root, const char *ctx)
and Is_Real_Container< NumContainer > and static Is_Real_Container< DenContainer > Array< Real > lfilter(const SignalContainer &signal, const NumContainer &numerator, const DenContainer &denominator, const Array< Real > &initial_state={})
static Array< Complex > prepare_stft_frame_input(const Array< Real > &frame, const Array< Real > &window, const size_t fft_size)
static Array< Array< Array< Complex > > > pbatched_stft(ThreadPool &pool, const Array< Array< Real > > &signals, const Array< Real > &window, const STFTOptions &options, const size_t chunk_size=0)
Parallel batched STFT across signals.
static Array< Array< Complex > > stft(const Container &signal, const size_t frame_size, const size_t hop_size, const bool pad_end=true)
static Array< Real > filtfilt(const Array< Real > &signal, const Array< Real > &coeffs, const size_t block_size=0)
Zero-phase FIR filtering via forward-backward convolution.
static Array< T > to_array(const Container &input)
static Array< Complex > compact_real_spectrum(const Array< Complex > &full_spectrum)
static Array< size_t > axis_base_offsets(const Array< size_t > &shape, const Array< size_t > &strides, const size_t axis)
Rationale: Generates the base memory offsets for all 1D slices along a specific axis in a multidimens...
static Real kaiser_beta(const Real attenuation_db)
Returns the Kaiser beta that corresponds to an attenuation goal.
static Array< PoleZeroPair > pair_poles_and_zeros(const BiquadSection §ion)
static Array< Complex > pspectrum(ThreadPool &pool, const Array< Real > &input, const size_t chunk_size=0)
Parallel version of spectrum(const Array<Real>&).
static Array< T > slice_copy(const Array< T > &input, const size_t offset, const size_t length)
Functional slice extraction.
static Array< PoleZeroPair > pair_poles_and_zeros(const Array< Complex > &zeros, const Array< Complex > &poles)
Greedily pairs zeros and poles by nearest distance.
static void transform_axes(Array< Complex > &data, const TensorLayout &layout, const Array< size_t > &axes, const bool invert)
In-place FFT/IFFT along multiple tensor axes of a flat buffer.
static Array< size_t > tensor3_shape(const Array< Array< Array< Complex > > > &input, const char *ctx)
Rationale: Validates that a nested Array represents a 3D tensor and returns its {d0,...
and static Is_Real_Container< CoeffContainer > Array< Real > filtfilt(const SignalContainer &signal, const CoeffContainer &coeffs, const size_t block_size=0)
and static Is_Biquad_Container< SectionsContainer > Array< Array< Real > > batched_sosfilt(const SignalsContainer &signals, const SectionsContainer §ions)
static Array< BiquadSection > transfer_function_to_sections(const Array< Real > &numerator, const Array< Real > &denominator, const char *ctx)
overload of transfer_function_to_sections using raw coefficients.
static Array< Array< Real > > batched_sosfilt(const Array< Array< Real > > &signals, const Array< BiquadSection > §ions)
One-shot batched causal SOS filtering across real channels.
static const char * simd_backend_name() noexcept
Returns the active default SIMD backend name for this precision.
static Array< Array< Real > > pirfft_batch(ThreadPool &pool, const Array< Array< Complex > > &spectra, const size_t signal_size, const size_t chunk_size=0)
Parallel compact real inverse FFT for equal-length batches.
static Array< Real > magnitude_spectrum(const Array< Complex > &input)
Returns |X[k]| for each frequency bin in an FFT output.
static Array< Complex > initialize_root_guesses(const Array< Real > &coefficients)
static Array< Real > phase_delay(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points=512, const bool whole=false)
static FrequencyResponse freqz(const SectionsContainer §ions, const size_t num_points=512, const bool whole=false)
static Array< Complex > polynomial_from_roots_complex(const Array< Complex > &roots)
Rationale: Polynomial expansion (z-r1)(z-r2)...
static Array< Real > multiply(const Array< Real > &a, const Array< Real > &b)
Multiplies two real-valued sequences using convolution.
and static Is_Complex_Container< Container2 > Array< Complex > multiply(const Container1 &a, const Container2 &b)
Multiplication for complex-valued containers.
static void validate_stable(const IIRCoefficients &coeffs)
and static Is_Real_Container< DenContainer > void validate_no_near_pole_zero_cancellation(const NumContainer &numerator, const DenContainer &denominator, const Real tolerance)
static AnalogPrototype chebyshev2_prototype(const size_t order, const Real attenuation_db, const char *ctx)
static Array< Complex > transform(const Array< Real > &input)
Computes the FFT of a real-valued input sequence.
and Is_Real_Container< NumContainer > and static Is_Real_Container< DenContainer > Array< Array< Real > > batched_lfilter(const SignalsContainer &signals, const NumContainer &numerator, const DenContainer &denominator, const Array< Array< Real > > &initial_states={})
static Array< Complex > ptransform_padded(ThreadPool &pool, const Array< Complex > &input, const size_t chunk_size=0)
static Array< RootGroup > root_groups_impl(const Array< Complex > &raw_roots, const char *ctx)
Rationale: Partitions a set of roots into real sections and conjugate pairs to facilitate SOS decompo...
static Array< Array< Complex > > project_to_plan_batch_spectra(const Plan &plan, const Array< Array< Real > > &input, ThreadPool *pool=nullptr, const size_t chunk_size=0)
static AnalogPrototype bessel_prototype(const size_t order, const char *ctx)
implementation of the Bessel filter analog prototype.
static size_t next_power_of_two(const size_t n)
Calculates the smallest power of two greater than or equal to n.
static Array< BiquadSection > elliptic_lowpass(const size_t order, const Real ripple_db, const Real attenuation_db, const Real cutoff_frequency, const Real sample_rate)
Digital elliptic/Cauer low-pass design returned as SOS.
static PhaseMarginInfo phase_margin(const Array< BiquadSection > §ions, const size_t num_points=1024, const bool whole=false)
static Array< Real > pistft(ThreadPool &pool, const Array< Array< Complex > > &spectrogram, const Array< Real > &window, const ISTFTOptions &options, const size_t chunk_size=0)
Parallel ISTFT using a shared analysis/synthesis window.
static Array< Array< Complex > > stft_impl(const Array< Real > &signal, const Array< Real > &window, const STFTOptions &options, ThreadPool *pool=nullptr, const size_t chunk_size=0)
static Array< Real > real_polynomial_from_roots(const Array< Complex > &roots, const char *ctx)
Rationale: Constructs a real-coefficient polynomial from a set of complex roots, verifying that roots...
static const char * batched_plan_simd_backend_name() noexcept
Returns the active batch-plan SIMD backend name.
static Array< Complex > ptransform_padded(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of transform_padded(const Container&).
static bool has_near_pole_zero_cancellation(const SectionsContainer §ions, const Real tolerance)
static Array< Real > kaiser_window(const size_t n, const Real beta)
Returns a Kaiser window of length n.
static Array< Real > group_delay(const IIRCoefficients &coeffs, const size_t num_points=512, const bool whole=false)
static Array< T > polynomial_multiply(const Array< T > &lhs, const Array< T > &rhs)
Rationale: Standard polynomial multiplication (O(N*M)).
static Array< Real > power_spectrum(const Array< Complex > &input)
Returns |X[k]|^2 for each frequency bin in an FFT output.
static Array< Array< Complex > > inverse_transform_batch(const Array< Array< Complex > > &input)
Functional batch IFFT wrapper.
static Array< size_t > factor_small_radices(size_t n)
Rationale: Factors a number into small primes (2, 3, 4, 5) to determine if a composite FFT kernel is ...
static void transform_axis_impl(Array< Complex > &data, const TensorLayout &layout, const size_t axis, const bool invert, ThreadPool *pool=nullptr, const size_t chunk_size=0)
Rationale: Internal implementation of a 1D FFT along a single tensor axis.
static Array< Complex > poles(const BiquadSection §ion)
static Array< Real > resample_poly(const Array< Real > &signal, const size_t up, const size_t down, const ResamplePolyOptions &options={})
Polyphase resampling with an internally designed Kaiser FIR.
static Array< Complex > poles(const Container &denominator)
static PolynomialEvaluation evaluate_polynomial_and_derivative(const Array< Real > &coefficients, const Complex &x, const char *ctx)
static Array< Real > iir_filtfilt_impl(const Array< Real > &signal, const Array< Real > &numerator, const Array< Real > &denominator, const char *ctx)
static Array< Complex > pspectrum(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of spectrum(const Container&).
static Array< PoleZeroPair > pair_poles_and_zeros(const IIRCoefficients &coeffs)
and static Is_Real_Container< DenContainer > Array< Real > group_delay(const NumContainer &numerator, const DenContainer &denominator, const size_t num_points=512, const bool whole=false)
static Real transform_stages(const size_t n) noexcept
Returns log2(N) stages, used for numerical error scaling.
static Array< PoleZeroPair > pole_zero_pairs_impl(const Array< Complex > &zeros, const Array< Complex > &poles)
Rationale: Matches zeros and poles into pairs by minimizing Euclidean distance, aiding in SOS decompo...
static Array< Real > firwin_lowpass(const size_t num_taps, const Real cutoff_frequency, const Real sample_rate, const Array< Real > &window)
FIR low-pass design via the window method.
static Real window_coherent_gain(const Array< Real > &window)
Returns the average window gain.
static void drop_prefix(Array< T > &input, const size_t count)
and static Is_Real_Container< DenContainer > Array< PoleZeroPair > pair_poles_and_zeros(const NumContainer &numerator, const DenContainer &denominator)
static Array< Real > istft(const Array< Array< Complex > > &spectrogram, const size_t frame_size, const size_t hop_size, const size_t signal_length=0)
Reconstructs a real signal using a Hann window pair.
static Array< Real > firwin_bandpass(const size_t num_taps, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate, const Array< Real > &window)
FIR band-pass design via the window method.
static size_t count_non_empty_real_batch(const Array< Array< Real > > &batch) noexcept
static Real solve_elliptic_selectivity_modulus(const Real k1, const size_t order, const char *ctx)
static bool should_parallelize_batch_work(ThreadPool *pool, const size_t batch_size, const size_t transform_size, const size_t min_work=8192) noexcept
static Real carlson_rf_impl(Real x, Real y, Real z)
static bool satisfies_nola(const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size)
Returns true when the window pair satisfies NOLA.
static std::pair< Array< Real >, Array< Complex > > extract_repeated_unit_roots(const Array< Real > &coefficients, const char *ctx)
static const char * simd_preference_name() noexcept
Returns the requested SIMD policy name.
static PowerSpectralDensity welch(const Array< Real > &signal, const Array< Real > &window, const Real sample_rate, const WelchOptions &options={})
Welch one-sided PSD estimate using a custom window.
static Array< Real > iir_steady_state(const Array< Real > &numerator, const Array< Real > &denominator, const char *ctx)
static Array< BiquadSection > chebyshev2_bandpass(const size_t order, const Real attenuation_db, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Chebyshev-II band-pass design returned as SOS.
static Array< Complex > apply_hamming_window(const Container &signal)
and static Is_Real_Container< WindowContainer > Array< Real > apply_window(const SignalContainer &signal, const WindowContainer &window)
static Real window_energy(const Array< Real > &window)
Returns the sum of squared window samples.
and static Is_Real_Container< Container2 > Array< Real > partitioned_convolution(const Container1 &signal, const Container2 &kernel, const size_t partition_size=0)
Partitioned convolution for real-valued containers.
static GainMarginInfo gain_margin(const Array< BiquadSection > §ions, const size_t num_points=1024, const bool whole=false)
static Complex elliptic_cd_minus_imaginary(const Real u, const Real v, const Real modulus, const char *ctx)
static Real scaled_tolerance(const Real reference, const Real multiplier=Real(256)) noexcept
static size_t tensor_element_count(const Array< size_t > &shape, const char *ctx)
static Array< Complex > spectrum(const Array< Real > &input)
DSP alias for the real forward FFT.
static PhaseMarginInfo phase_margin(const IIRCoefficients &coeffs, const size_t num_points=1024, const bool whole=false)
and Is_Real_Container< DesiredContainer > and static Is_Real_Container< WeightContainer > Array< Real > firls(const size_t num_taps, const BandContainer &bands, const DesiredContainer &desired, const Real sample_rate, const WeightContainer &weights)
static Real minimum_pole_zero_distance(const IIRCoefficients &coeffs)
static GainMarginInfo gain_margin_impl(const FrequencyResponse &response)
implementation helper for Gain Margin calculation.
static Array< Complex > flatten_tensor3_row_major(const Array< Array< Array< Complex > > > &input, const char *ctx)
Rationale: Flattens a 3D tensor into a 1D row-major array.
static bool is_one_sided_interior_bin(const size_t bin, const size_t fft_size) noexcept
Identifies non-redundant spectrum bins (excluding DC and Nyquist).
and static Is_Complex_Container< Container2 > Array< Complex > pmultiply(ThreadPool &pool, const Container1 &a, const Container2 &b, const size_t chunk_size=0)
Parallel version of multiply(const Container1&, const Container2&).
static Array< T > reverse_copy(const Array< T > &input)
Functional reverse.
static Array< Complex > transform(const Container &input)
Forward FFT for real-valued containers.
static void enforce_real_polynomial_symmetry(Array< Complex > &roots, const Real tol) noexcept
static Array< Real > group_delay_impl(const FrequencyResponse &response)
implementation of numeric group delay calculation.
static Real minimum_pole_zero_distance(const SectionsContainer §ions)
static Array< Real > blackman_window(const size_t n)
Returns a Blackman window of length n.
static size_t default_filtfilt_pad_length(const size_t signal_size, const size_t coeff_size) noexcept
Rationale: Rule-of-thumb for filtfilt padding to allow transients to decay.
and static Is_Real_Container< DenContainer > bool has_near_pole_zero_cancellation(const NumContainer &numerator, const DenContainer &denominator, const Real tolerance)
static Array< Complex > inverse_transform(const Array< Complex > &input)
Computes the Inverse Fast Fourier Transform (IFFT).
static void polish_roots_with_newton(const Array< Real > &coefficients, Array< Complex > &roots, const Real tol, const size_t iterations) noexcept
static Array< Complex > spectrum(const Container &input)
DSP alias for real-valued container forward FFT.
static AnalogPrototype elliptic_prototype(const size_t order, const Real ripple_db, const Real attenuation_db, const char *ctx)
static Real sum_squares(const Array< Real > &input) noexcept
L2 norm squared.
static void ptransform_batch(ThreadPool &pool, Array< Array< Complex > > &batch, const bool invert, const size_t chunk_size=0)
Parallel in-place batch transform for equal-length inputs.
static Array< Real > zero_padded_copy(const Array< Real > &input, const size_t n)
Returns a copy of input zero-padded to size n.
static void validate_no_near_pole_zero_cancellation(const Array< Complex > &zeros, const Array< Complex > &poles, const Real tolerance)
static PhaseMarginInfo phase_margin(const FrequencyResponse &response)
Estimates phase margin around the unity-gain crossover.
static Array< Complex > transform_real_optimized(const Array< Real > &input, ThreadPool *pool=nullptr, const size_t chunk_size=0)
Rationale: Optimized FFT for real signals of size N (power-of-two).
and static Is_Real_Container< Container2 > Array< Real > overlap_add_convolution(const Container1 &signal, const Container2 &kernel, const size_t block_size=0)
Overlap-add convolution for real-valued containers.
static Real minimum_pole_zero_distance(const Array< Real > &numerator, const Array< Real > &denominator)
static Array< BiquadSection > chebyshev1_bandpass(const size_t order, const Real ripple_db, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Chebyshev-I band-pass design returned as SOS.
static Real series_zero_tolerance(const size_t coeff_count, const Real multiplier=Real(256)) noexcept
static Array< Complex > transformed(const Array< Complex > &input, const bool invert=false)
Computes the FFT or IFFT and returns a new array.
static constexpr bool Is_Real_Container
static bool overlap_profile_has_nola(const Array< Real > &profile) noexcept
Rationale: Non-zero Overlap-Add (NOLA) ensures the signal can be reconstructed (it's never multiplied...
static Real integrate_linear_cos_basis(const Real omega_lo, const Real omega_hi, const Real desired_lo, const Real desired_hi, const size_t harmonic) noexcept
implementation helper for firls linear band integration.
static bool has_near_pole_zero_cancellation(const IIRCoefficients &coeffs, const Real tolerance)
static IIRCoefficients normalize_iir_coefficients(const Array< Real > &numerator, const Array< Real > &denominator, const char *ctx)
static Array< Array< Complex > > transformed_batch(const Array< Array< Complex > > &input, const bool invert=false)
Functional batch FFT/IFFT wrapper.
static Array< Real > istft(const Array< Array< Complex > > &spectrogram, const size_t frame_size, const ISTFTOptions &options)
Reconstructs a Hann-window STFT with explicit options.
static bool try_laguerre_roots(const Array< Real > &coefficients, const Real tol, Array< Complex > &roots) noexcept
static Array< BiquadSection > design_prototype_sections(const AnalogPrototype &prototype, const size_t order, const Real cutoff_frequency, const Real sample_rate, const bool highpass, const char *ctx)
static FrequencyResponse freqz(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points=512, const bool whole=false)
Samples the discrete-time transfer response on the unit circle.
static Real max_abs_value(const Array< Real > &input) noexcept
L-infinity norm.
static Array< Real > bilinear_substitute_polynomial(const Array< Real > &analog, const size_t order, const Real sample_rate, const char *ctx)
Rationale: Performs the bilinear substitution mapping from the S-plane to the Z-domain: H(z) = H(s)|s...
static Array< Complex > zero_padded_copy(const Array< Complex > &input, const size_t n)
Returns a copy of input zero-padded to size n.
static bool has_near_pole_zero_cancellation(const BiquadSection §ion, const Real tolerance)
static Real integrate_cos_basis(const Real omega_lo, const Real omega_hi, const size_t harmonic) noexcept
Rationale: Evaluates integral of cos(k*w) from omega_lo to omega_hi.
static void validate_no_near_pole_zero_cancellation(const Array< Real > &numerator, const Array< Real > &denominator, const Real tolerance)
static Array< BiquadSection > design_low_high_sections(const AnalogPrototype &prototype, const Real cutoff_frequency, const Real sample_rate, const bool highpass, const char *ctx)
internal helper for designing standard lowpass/highpass SOS filters.
static Array< BiquadSection > chebyshev1_highpass(const size_t order, const Real ripple_db, const Real cutoff_frequency, const Real sample_rate)
Digital Chebyshev-I high-pass design returned as SOS.
and static Is_Real_Container< DesiredContainer > Array< Real > firls(const size_t num_taps, const BandContainer &bands, const DesiredContainer &desired, const Real sample_rate)
static CrossSpectralDensity csd(const Array< Real > &x, const Array< Real > &y, const size_t frame_size, const Real sample_rate, const WelchOptions &options={})
One-sided cross-spectral density estimate using a Hann window.
static Array< Real > firwin_lowpass_impl(const size_t num_taps, const Real cutoff_frequency, const Real sample_rate, const Array< Real > &window, const char *ctx)
implementation of the window-method FIR design.
and static Is_Real_Container< WindowContainer > Array< Complex > apply_window(const SignalContainer &signal, const WindowContainer &window)
static void apply_one_sided_density_scaling(Density &value, const size_t bin, const size_t fft_size) noexcept
Rationale: Multiplies spectral densities by 2 for interior bins to account for energy in the omitted ...
static Array< Real > pmultiply(ThreadPool &pool, const Array< Real > &a, const Array< Real > &b, const size_t chunk_size=0)
Parallel version of multiply(const Array<Real>&, const Array<Real>&).
static Array< Complex > ptransform_padded(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel zero-pad and forward-FFT for a generic complex iterable.
static Array< Array< Array< Complex > > > multichannel_stft(const Array< Array< Real > > &signals, const Array< Real > &window, const STFTOptions &options={}, const SpectrogramLayout layout=SpectrogramLayout::channel_frame_bin)
Multichannel STFT with explicit channel-major or frame-major output layout.
static Real minimum_pole_zero_distance(const BiquadSection §ion)
static bool satisfies_cola(const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size)
Returns true when the window pair satisfies COLA.
static void validate_no_near_pole_zero_cancellation(const SectionsContainer §ions, const Real tolerance)
and static Is_Real_Container< DenContainer > FrequencyResponse freqz(const NumContainer &numerator, const DenContainer &denominator, const size_t num_points=512, const bool whole=false)
static bool neon_runtime_available() noexcept
Returns whether NEON dispatch is supported at runtime.
and static Is_Real_Container< WindowContainer > Array< Complex > windowed_spectrum(const SignalContainer &signal, const WindowContainer &window)
static Array< Real > substitute_rational_polynomial(const Array< Real > &poly, const Array< Real > &numerator, const Array< Real > &denominator, const char *ctx)
Rationale: Evaluates P(N(z)/D(z)) by scaling by the common denominator D(z)^degree and evaluating the...
static TransferTerms evaluate_transfer_terms_at(const Array< Real > &numerator, const Array< Real > &denominator, const Real omega, const char *ctx)
implementation helper for evaluating H(z) and its derivatives on the unit circle.
static AnalogPrototype chebyshev1_prototype(const size_t order, const Real ripple_db, const char *ctx)
static Array< Complex > ptransform(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of transform(const Container&).
static Array< Real > iir_filter_impl(const Array< Real > &signal, const Array< Real > &numerator, const Array< Real > &denominator, const Array< Real > &initial_state, const char *ctx, Array< Real > *final_state=nullptr)
static void validate_stable(const BiquadSection §ion)
static Array< Complex > windowed_spectrum(const Array< Real > &signal, const Array< Real > &window)
Returns the FFT of a real signal after applying a window.
static Array< size_t > choose_remez_extrema(const Array< Real > &weighted_error, const size_t required_count, const char *ctx)
static Array< Complex > transformed_axes(const Array< Complex > &input, const TensorLayout &layout, const Array< size_t > &axes, const bool invert=false)
Functional tensor FFT/IFFT wrapper for flat buffers.
static Array< Complex > ptransform_padded(ThreadPool &pool, const Array< Real > &input, const size_t chunk_size=0)
Parallel version of transform_padded(const Array<Real>&).
and static Is_Real_Container< DenContainer > PhaseMarginInfo phase_margin(const NumContainer &numerator, const DenContainer &denominator, const size_t num_points=1024, const bool whole=false)
static GainMarginInfo gain_margin(const IIRCoefficients &coeffs, const size_t num_points=1024, const bool whole=false)
static void validate_istft_configuration(const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t fft_size, const ISTFTOptions &options, const char *ctx)
Rationale: Validates and normalizes ISTFT parameters, checking window constraints and FFT size.
static Array< Real > firwin_highpass(const size_t num_taps, const Real cutoff_frequency, const Real sample_rate, const Array< Real > &window)
FIR high-pass design via spectral inversion of a low-pass design.
static GainMarginInfo gain_margin(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points=1024, const bool whole=false)
and static Is_Real_Container< DenContainer > IIRCoefficients bilinear_transform(const NumContainer &analog_numerator, const DenContainer &analog_denominator, const Real sample_rate)
static Array< Real > group_delay(const Container &numerator, const size_t num_points=512, const bool whole=false)
static Array< Complex > enforce_conjugate_symmetry(const Array< Complex > &roots)
Enforce conjugate symmetry on roots of a real-coefficient polynomial.
static Real real_projection_tolerance(const Complex &value, const size_t n) noexcept
Rationale: Estimates the maximum expected numerical noise for a value after an N-point transform,...
static Array< Complex > prfft(ThreadPool &pool, const Array< Real > &input, const size_t chunk_size=0)
Parallel version of rfft(const Array<Real>&).
static Array< Array< Complex > > reshape_matrix_row_major(const Array< Complex > &input, const size_t rows, const size_t cols, const char *ctx)
Rationale: Reshapes a 1D row-major array back into a 2D matrix.
static void ptransform_axis(ThreadPool &pool, Array< Complex > &data, const TensorLayout &layout, const size_t axis, const bool invert, const size_t chunk_size=0)
Parallel in-place 1-D FFT/IFFT along one tensor axis.
static Array< Complex > poles(const SectionsContainer §ions)
and Is_Real_Container< NumContainer > and static Is_Real_Container< DenContainer > Array< Real > filtfilt(const SignalContainer &signal, const NumContainer &numerator, const DenContainer &denominator)
static void transform_impl(Array< Complex > &a, const bool invert, ThreadPool *pool=nullptr, const size_t chunk_size=0)
implementation of the power-of-two iterative FFT using mixed radix-2/4.
static Array< BiquadSection > chebyshev2_bandstop(const size_t order, const Real attenuation_db, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Chebyshev-II band-stop design returned as SOS.
static bool is_stable(const BiquadSection §ion)
static Array< Real > pinverse_transform_real(ThreadPool &pool, const Array< Complex > &input, const size_t chunk_size=0)
Parallel version of inverse_transform_real(const Array<Complex>&).
static Array< Complex > ptransform(ThreadPool &pool, const Array< Real > &input, const size_t chunk_size=0)
Parallel real-input FFT.
static Real sum_values(const Array< Real > &input) noexcept
Sum reduction.
static FrequencyResponse freqz(const Array< Real > &numerator, const size_t num_points=512, const bool whole=false)
FIR frequency response.
static Array< Real > cosine_sum_window(const size_t n, const Real a0, const Real a1=Real(0), const Real a2=Real(0))
Internal generator for generalized cosine windows (Hann, Hamming, Blackman).
static void validate_stable(const Array< BiquadSection > §ions, const Real min_margin)
static void transform_batch(Array< Array< Complex > > &batch, const bool invert)
In-place batch transform for equal-length complex inputs.
static void transform_axis(Array< Complex > &data, const TensorLayout &layout, const size_t axis, const bool invert)
In-place 1-D FFT/IFFT along one tensor axis of a flat buffer.
static Array< Array< Array< Complex > > > pmultichannel_stft(ThreadPool &pool, const Array< Array< Real > > &signals, const Array< Real > &window, const STFTOptions &options={}, const SpectrogramLayout layout=SpectrogramLayout::channel_frame_bin, const size_t chunk_size=0)
Parallel multichannel STFT with selectable output layout.
static bool satisfies_cola(const Array< Real > &window, const size_t hop_size)
Returns true when one window satisfies COLA by itself.
static GainMarginInfo gain_margin_refined_impl(const FrequencyResponse &response, const Evaluator &evaluator)
implementation of refined Gain Margin calculation using analytic evaluation.
static bool has_near_pole_zero_cancellation(const Array< BiquadSection > §ions, const Real tolerance)
static void validate_stable(const SectionsContainer §ions, const Real min_margin)
static void validate_stable(const Array< Real > &denominator)
static Array< Complex > ptransform(ThreadPool &pool, const Container &input, const size_t chunk_size=0)
Parallel version of transform(const Container&).
static Array< Real > apply_window(const Array< Real > &signal, const Array< Real > &window)
Applies a real window sample-by-sample to a real signal.
static Real complete_elliptic_first_kind(const Real modulus, const char *ctx)
Rationale: Complete elliptic integral K(k) using the modulus convention required by elliptic filter d...
static PowerSpectralDensity welch(const Array< Real > &signal, const size_t frame_size, const Real sample_rate, const WelchOptions &options={})
Welch one-sided PSD estimate using a Hann window.
static Real integer_power(const Real base, const size_t exponent) noexcept
Fast integer power utility.
static Array< Real > group_delay(const Array< Real > &numerator, const Array< Real > &denominator, const size_t num_points=512, const bool whole=false)
static Array< BiquadSection > chebyshev1_bandstop(const size_t order, const Real ripple_db, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Chebyshev-I band-stop design returned as SOS.
static Array< Array< Real > > batched_istft(const Array< Array< Array< Complex > > > &spectrograms, const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const ISTFTOptions &options, const Array< size_t > &signal_lengths={})
Batched ISTFT with optional per-signal output lengths.
static size_t tensor_max_offset(const Array< size_t > &shape, const Array< size_t > &strides, const char *ctx)
static Array< Real > poverlap_add_convolution(ThreadPool &pool, const Array< Real > &signal, const Array< Real > &kernel, const size_t block_size=0, const size_t chunk_size=0)
Parallel convenience wrapper for overlap-add convolution.
and static Is_Real_Container< DenContainer > Array< Real > phase_delay(const NumContainer &numerator, const DenContainer &denominator, const size_t num_points=512, const bool whole=false)
static PhaseMarginInfo phase_margin(const SectionsContainer §ions, const size_t num_points=1024, const bool whole=false)
static Array< Real > remez_impl(const size_t num_taps, const Array< Real > &bands, const Array< Real > &desired, const Real sample_rate, const Array< Real > &weights, const size_t grid_density, const size_t max_iterations, const char *ctx)
static Array< Complex > ptransformed(ThreadPool &pool, const Container &input, const bool invert=false, const size_t chunk_size=0)
Parallel FFT/IFFT for a generic complex iterable.
static bool is_stable(const IIRCoefficients &coeffs)
static Array< BiquadSection > bessel_bandpass(const size_t order, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Bessel band-pass design returned as SOS.
static void validate_stable(const Container &denominator, const Real min_margin)
static Array< Real > apply_hann_window(const Container &signal)
static IIRCoefficients bilinear_transform_impl(const Array< Real > &analog_numerator, const Array< Real > &analog_denominator, const Real sample_rate, const char *ctx)
implementation helper for the bilinear transform.
static Array< BiquadSection > design_bandstop_sections_without_numerator_roots(const AnalogPrototype &prototype, const size_t order, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate, const char *ctx)
Specialization for band-stop designs to handle zeros explicitly.
and static Is_Real_Container< WindowContainer > Array< Array< Complex > > pstft(ThreadPool &pool, const SignalContainer &signal, const WindowContainer &window, const STFTOptions &options, const size_t chunk_size=0)
static Array< size_t > row_major_strides(const Array< size_t > &shape, const char *ctx)
SpectrogramLayout
Logical layout for multichannel spectrograms.
static Array< Real > phase_delay(const Container &numerator, const size_t num_points=512, const bool whole=false)
static Array< Complex > transformed(const Container &input, const bool invert=false)
FFT/IFFT for a generic complex iterable.
static Array< Complex > transform_padded(const Array< Real > &input)
Computes the FFT after zero-padding the input to the next power of two.
static void transform(Array< Complex > &a, const bool invert)
Computes the Fast Fourier Transform (FFT) in-place.
static Real window_enbw(const Array< Real > &window)
Returns the equivalent noise bandwidth of a window in bins.
static Array< Array< Complex > > ptransformed_batch(ThreadPool &pool, const Array< Array< Complex > > &input, const bool invert=false, const size_t chunk_size=0)
Parallel functional batch FFT/IFFT wrapper.
static Array< Real > phase_delay(const IIRCoefficients &coeffs, const size_t num_points=512, const bool whole=false)
static Array< Complex > project_to_plan_spectrum(const Plan &plan, const Array< Real > &input, ThreadPool *pool=nullptr, const size_t chunk_size=0)
static Real prewarp_frequency(const Real cutoff_frequency, const Real sample_rate, const char *ctx)
Rationale: Prewarps an analog cutoff frequency to compensate for the non-linear frequency mapping of ...
static Array< Real > apply_hamming_window(const Array< Real > &signal)
Applies a Hamming window of matching size to a real signal.
static JacobiValues jacobi_sn_cn_dn(const Real argument, const Real modulus, const char *ctx)
Rationale: Simultanous calculation of Jacobi sn, cn, dn functions using the descending Landen transfo...
static Array< size_t > matrix_shape(const Array< Array< Complex > > &input, const char *ctx)
Rationale: Validates that an Array of Arrays represents a rectangular matrix and returns its {rows,...
static Array< Complex > flatten_matrix_row_major(const Array< Array< Complex > > &input, const char *ctx)
Rationale: Flattens a 2D matrix into a 1D row-major array.
static Array< Real > partitioned_convolution(const Array< Real > &signal, const Array< Real > &kernel, const size_t partition_size=0)
Convenience wrapper for low-latency partitioned convolution.
static Array< Array< Array< Complex > > > transpose_spectrogram_layout(const Array< Array< Array< Complex > > > &input, const SpectrogramLayout source, const SpectrogramLayout target)
Converts multichannel spectrograms between channel-major and frame-major layouts.
static Array< Real > filtfilt(const Array< Real > &signal, const Array< Real > &numerator, const Array< Real > &denominator)
Zero-phase IIR filtering using transfer-function coefficients.
static Array< Complex > expand_real_spectrum(const Array< Complex > &spectrum, const size_t signal_size, const char *ctx)
static void transform_any_size_impl(Array< Complex > &a, const bool invert, ThreadPool *pool=nullptr, const size_t chunk_size=0)
static Array< Array< Complex > > transformed2d(const Array< Array< Complex > > &input, const bool invert=false)
Functional 2-D FFT/IFFT wrapper for rectangular complex matrices.
static Array< Real > pistft(ThreadPool &pool, const Array< Array< Complex > > &spectrogram, const Array< Real > &window, const size_t hop_size, const size_t signal_length=0, const size_t chunk_size=0)
Parallel STFT inversion using a shared analysis/synthesis window.
static Array< Real > inverse_transform_real(const Array< Complex > &input)
Computes the IFFT and projects the result back to real values.
static Array< Complex > pmultiply(ThreadPool &pool, const Array< Complex > &a, const Array< Complex > &b, const size_t chunk_size=0)
Parallel version of multiply(const Array<Complex>&, const Array<Complex>&).
static constexpr bool Is_Complex_Container
static Array< Real > scaled_copy(const Array< Real > &input, const Real factor)
Scalar scaling.
static Array< Complex > gather_axis_slice(const Array< Complex > &data, const size_t base_offset, const size_t axis_length, const size_t axis_stride)
Rationale: Extracts a non-contiguous slice of data along a tensor axis into a contiguous array for FF...
static Array< Real > istft(const Array< Array< Complex > > &spectrogram, const Array< Real > &window, const size_t hop_size, const size_t signal_length=0)
Reconstructs a real signal using the same window for analysis and synthesis.
static Array< Array< Array< Complex > > > transformed3d(const Array< Array< Array< Complex > > > &input, const bool invert=false)
Functional 3-D FFT/IFFT wrapper for rectangular complex tensors.
static size_t recommended_cache_tile_size(const size_t transform_size, const size_t batch_size, ThreadPool *pool) noexcept
static Array< Complex > pinverse_transform(ThreadPool &pool, const Array< Complex > &input, const size_t chunk_size=0)
Parallel version of inverse_transform(const Array<Complex>&).
static IIRCoefficients bilinear_transform(const Array< Real > &analog_numerator, const Array< Real > &analog_denominator, const Real sample_rate)
Bilinear transform of an analog transfer function.
static size_t resolve_welch_hop_size(const WelchOptions &options, const size_t frame_size, const char *ctx)
Default hop size for Welch analysis (50% overlap).
and static Is_Real_Container< CoeffContainer > Array< Real > upfirdn(const SignalContainer &signal, const CoeffContainer &coeffs, const size_t up=1, const size_t down=1)
static SimdPreference simd_preference() noexcept
Returns the runtime SIMD preference requested via environment.
static Array< Complex > apply_hann_window(const Container &signal)
static void trim_to_size(Array< Complex > &input, const size_t n)
Trims input to size n.
static void validate_stable(const BiquadSection §ion, const Real min_margin)
static Real polynomial_root_lower_bound(const Array< Real > &monic) noexcept
and static Is_Biquad_Container< SectionsContainer > Array< Real > filtfilt(const SignalContainer &signal, const SectionsContainer §ions)
static Array< Array< Array< Complex > > > reshape_tensor3_row_major(const Array< Complex > &input, const size_t dim0, const size_t dim1, const size_t dim2, const char *ctx)
Rationale: Reshapes a 1D row-major array back into a 3D tensor.
static Real stability_margin(const BiquadSection §ion)
static Array< BiquadSection > butterworth_bandstop(const size_t order, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate)
Digital Butterworth band-stop design returned as SOS.
static PhaseMarginInfo phase_margin_refined_impl(const FrequencyResponse &response, const Evaluator &evaluator)
implementation of refined Phase Margin calculation using analytic evaluation.
static Array< Real > firls(const size_t num_taps, const Array< Real > &bands, const Array< Real > &desired, const Real sample_rate, const Array< Real > &weights={})
FIR design by weighted least squares over piecewise-linear bands.
static Array< Array< Complex > > inverse_transform2d(const Array< Array< Complex > > &input)
Functional inverse 2-D FFT wrapper.
static Array< Array< Real > > pbatched_istft(ThreadPool &pool, const Array< Array< Array< Complex > > > &spectrograms, const Array< Real > &window, const ISTFTOptions &options, const Array< size_t > &signal_lengths={}, const size_t chunk_size=0)
Parallel batched ISTFT using a shared analysis/synthesis window.
static Array< Array< Real > > irfft_batch(const Array< Array< Complex > > &spectra, const size_t signal_size)
Functional compact real inverse FFT for equal-length batches.
static Array< size_t > frame_offsets_impl(const size_t signal_size, const size_t frame_size, const size_t hop_size, const bool pad_end, const char *ctx)
Computes framing offsets.
and static Is_Biquad_Container< SectionsContainer > Array< Real > sosfilt(const SignalContainer &signal, const SectionsContainer §ions)
static Array< Real > zero_pad_edges(const Array< Real > &signal, const size_t left_pad, const size_t right_pad)
functional zero-padding at both ends.
static Real polynomial_root_upper_bound(const Array< Real > &monic) noexcept
static Array< PoleZeroPair > pair_poles_and_zeros(const SectionsContainer §ions)
static Array< Array< Complex > > pinverse_transform_batch(ThreadPool &pool, const Array< Array< Complex > > &input, const size_t chunk_size=0)
Parallel functional batch IFFT wrapper.
static Real integrate_cos_product(const Real omega_lo, const Real omega_hi, const size_t lhs_harmonic, const size_t rhs_harmonic) noexcept
Rationale: Evaluates integral of cos(m*w)*cos(n*w).
static Array< Real > reverse_bessel_polynomial(const size_t order, const char *ctx)
Rationale: Generates coefficients for the Reverse Bessel Polynomial used in Bessel filter prototypes.
static Array< Real > phase_spectrum(const Array< Complex > &input)
Returns arg(X[k]) for each frequency bin in an FFT output.
static Array< Real > magnitude_spectrum(const Container &input)
Magnitude spectrum for complex-valued containers.
static void ptransform_axes(ThreadPool &pool, Array< Complex > &data, const TensorLayout &layout, const Array< size_t > &axes, const bool invert, const size_t chunk_size=0)
Parallel in-place FFT/IFFT along multiple tensor axes.
static Array< size_t > normalize_axes(const Array< size_t > &axes, const size_t rank, const char *ctx)
Rationale: Validates axis indices against the tensor rank and ensures no axis is specified more than ...
static Array< Real > window_overlap_profile(const Array< Real > &analysis_window, const Array< Real > &synthesis_window, const size_t hop_size)
Returns the overlap-add normalization profile for one hop period.
static TensorLayout normalize_tensor_layout(const Array< Complex > &data, const TensorLayout &layout, const char *ctx)
Rationale: Validates a tensor layout and infers row-major strides if they are not provided.
static Array< Array< Complex > > transform_stft_frames(const Array< Array< Real > > &frames, const Array< Real > &window, const size_t fft_size, const Plan &plan, ThreadPool *pool, const size_t chunk_size)
static Array< Complex > transform_padded(const Container &input)
Zero-pad and forward-FFT a generic complex iterable.
static void validate_stable(const Container &denominator)
static Array< BiquadSection > chebyshev1_lowpass(const size_t order, const Real ripple_db, const Real cutoff_frequency, const Real sample_rate)
Digital Chebyshev-I low-pass design returned as SOS.
static Array< Array< Complex > > stft(const Array< Real > &signal, const size_t frame_size, const STFTOptions &options)
Computes a Hann-window STFT with explicit analysis options.
static Array< Array< Array< Complex > > > batched_stft(const Array< Array< Real > > &signals, const Array< Real > &window, const STFTOptions &options)
Batched STFT over a collection of real signals.
static Array< Real > filtfilt(const SignalContainer &signal, const IIRCoefficients &coeffs)
static Real inverse_jacobi_sc(const Real value, const Real modulus, const char *ctx)
static Array< Real > inverse_transform_real(const Container &input)
Inverse real transform for complex-valued containers.
static Array< Array< Array< Complex > > > transformed2d_batch(const Array< Array< Array< Complex > > > &input, const bool invert=false)
Functional 2-D batched FFT/IFFT wrapper over a matrix stack.
static Array< Complex > zeros(const Array< BiquadSection > §ions)
static Array< Real > group_delay(const BiquadSection §ion, const size_t num_points=512, const bool whole=false)
static Complex twiddle_at(const Real angle, const size_t index)
Computes a single twiddle factor exp(j * angle * index).
static Array< Real > reflect_pad_signal(const Array< Real > &signal, const size_t pad_len)
Rationale: Applies reflection padding at signal edges to reduce boundary artifacts during filtering.
std::complex< Real > Complex
static Array< Complex > apply_blackman_window(const Container &signal)
static Real evaluate_analog_transfer_magnitude(const Array< Real > &numerator, const Array< Real > &denominator, const Real omega, const char *ctx)
Rationale: Direct evaluation of |H(j*omega)| for an analog transfer function.
static Array< Real > overlap_add_frames(const Array< Array< Real > > &frames, const size_t hop_size, const size_t signal_length=0)
Overlap-adds a frame sequence with a fixed hop size.
static Array< Complex > transform_padded(const Array< Complex > &input)
static Array< Real > overlap_add_convolution(const Array< Real > &signal, const Array< Real > &kernel, const size_t block_size=0)
Convenience wrapper for long real convolution via overlap-add.
static Array< Complex > poles(const Array< Real > &denominator)
Returns the poles of a transfer denominator.
static bool satisfies_nola(const Array< Real > &window, const size_t hop_size)
Returns true when one window satisfies NOLA by itself.
static Array< Real > pistft(ThreadPool &pool, const Array< Array< Complex > > &spectrogram, const size_t frame_size, const ISTFTOptions &options, const size_t chunk_size=0)
Parallel Hann-window ISTFT with explicit options.
static bool is_stable(const Array< BiquadSection > §ions)
static bool overlap_profile_has_cola(const Array< Real > &profile) noexcept
Rationale: Constant Overlap-Add (COLA) ensures that OLA reconstruction has no amplitude modulation ar...
static Array< Real > firwin_bandstop(const size_t num_taps, const Real low_cutoff_frequency, const Real high_cutoff_frequency, const Real sample_rate, const Real attenuation_db)
FIR band-stop design using a Kaiser window.
static Array< Real > resample_poly(const Array< Real > &signal, const size_t up, const size_t down, const Array< Real > &coeffs)
Polyphase resampling with explicit FIR coefficients.
static void add_scaled_polynomial(Array< T > &dst, const Array< T > &src, const T scale)
Rationale: In-place addition of a scaled polynomial.
static Array< Complex > windowed_spectrum(const Array< Complex > &signal, const Array< Real > &window)
Returns the FFT of a complex signal after applying a window.
static size_t effective_coeff_length(const Array< Real > &input) noexcept
Returns length excluding trailing insignificant zeros.
static Array< Complex > multiply_complex_impl(const Array< Complex > &a, const Array< Complex > &b, ThreadPool *pool=nullptr, const size_t chunk_size=0)
Rationale: Linear convolution of two complex sequences via FFT.
Graph implemented with double-linked adjacency lists.
A reusable thread pool for efficient parallel task execution.
size_t num_threads() const noexcept
Get the number of worker threads.
Minimal std::expected-style result type for C++20.
__gmp_expr< T, __gmp_unary_expr< __gmp_expr< T, U >, __gmp_y1_function > > y1(const __gmp_expr< T, U > &expr)
__gmp_expr< T, __gmp_unary_expr< __gmp_expr< T, U >, __gmp_floor_function > > floor(const __gmp_expr< T, U > &expr)
__gmp_expr< typename __gmp_resolve_expr< T, V >::value_type, __gmp_binary_expr< __gmp_expr< T, U >, __gmp_expr< V, W >, __gmp_remainder_function > > remainder(const __gmp_expr< T, U > &expr1, const __gmp_expr< V, W > &expr2)
__gmp_expr< T, __gmp_unary_expr< __gmp_expr< T, U >, __gmp_y0_function > > y0(const __gmp_expr< T, U > &expr)
__gmp_expr< T, __gmp_binary_expr< __gmp_expr< T, U >, unsigned long int, __gmp_root_function > > root(const __gmp_expr< T, U > &expr, unsigned long int l)
size_t blossom_maximum_cardinality_matching(const GT &g, DynDlist< typename GT::Arc * > &matching, SA sa=SA())
Alias of compute_maximum_cardinality_general_matching().
const long double offset[]
Offset values indexed by symbol string length (bounded by MAX_OFFSET_INDEX)
constexpr State invert(const State &s) noexcept
Toggle a cell value between its zero and a non-zero counterpart.
double density(const Lattice &lat, const typename Lattice::state_type &s)
Compiler_SSA_Block & block(Compiler_SSA_Function &function, const Compiler_SSA_Block_Id id)
Main namespace for Aleph-w library functions.
size_t size(Node *root) noexcept
std::pair< TgtContainer< typename SrcContainer::Item_Type >, TgtContainer< typename SrcContainer::Item_Type > > partition(const SrcContainer &c, std::function< bool(const typename SrcContainer::Item_Type &)> operation)
Partition a container into two based on a predicate.
static long & low(typename GT::Node *p)
Internal helper: low-link value stored directly in NODE_COOKIE(p).
DynList< T > repeated(const Container< T > &c)
Return elements that appear more than once in the container.
void parallel_for_index(ThreadPool &pool, size_t start, size_t end, F &&f, size_t chunk_size=0)
Apply a function to each element in parallel (index-based).
and
Check uniqueness with explicit hash + equality functors.
std::decay_t< typename HeadC::Item_Type > T
auto mean(const Container &data) -> std::decay_t< decltype(*std::begin(data))>
Compute the arithmetic mean.
bool diff(const C1 &c1, const C2 &c2, Eq e=Eq())
Check if two containers differ.
std::pair< First, Second > pair
Alias to std::pair kept for backwards compatibility.
void error(const char *file, int line, const char *format,...)
Print an error message with file and line info.
void next()
Advance all underlying iterators (bounds-checked).
auto mode(const Container &data) -> std::decay_t< decltype(*std::begin(data))>
Compute the mode (most frequent value).
auto max_value(const Container &data) -> std::decay_t< decltype(*std::begin(data))>
Compute maximum value.
Itor::difference_type count(const Itor &beg, const Itor &end, const T &value)
Count elements equal to a value.
T sum(const Container &container, const T &init=T{})
Compute sum of all elements.
static struct argp_option options[]
static void section(const string &title)
Parameters for an analog filter prototype (S-domain).
A stable Second-Order Section (SOS) building block.
Array< Real > denominator() const
Array< Real > numerator() const
Magnitude-squared coherence sampled in Hertz.
Array< Real > magnitude_squared
One-sided cross-spectral density estimate sampled in Hertz.
Discrete frequency response sampled on a fixed angular grid.
Array< Complex > response
Array< Real > power() const
Array< Real > magnitude() const
Array< Real > phase() const
Gain-margin estimate around a phase crossover.
Coefficients for an Infinite Impulse Response (IIR) filter.
Array< Real > denominator
Feed-backward coefficients (a).
Array< Real > numerator
Feed-forward coefficients (b).
Options for ISTFT reconstruction.
Value triplet for Jacobi elliptic functions (sn, cn, dn).
Phase-margin estimate around a gain crossover.
Greedy nearest-neighbor pole/zero pairing entry.
bool is_cancellation(const Real tolerance) const noexcept
One-sided power spectral density estimate sampled in Hertz.
Options for default polyphase resampling filter design.
Rationale: Group of one or more roots (typically a conjugate pair) sharing a common geometric center.
Options for STFT analysis.
Layout descriptor for a flat tensor buffer.
Internal storage for complex transfer function evaluation terms.
Complex numerator_derivative
Complex denominator_derivative
Frequency-domain grid with target values and importance weights.
Options for Welch PSD, CSD, and coherence estimation.
FooMap m(5, fst_unit_pair_hash, snd_unit_pair_hash)
A modern, efficient thread pool for parallel task execution.
Dynamic array container with automatic resizing.