WindowFunctions.h

Go to the documentation of this file.

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <numbers>
 
namespace dspark {
 
template <typename T>
struct WindowFunctions
{
    // Prevent instantiation (static-only utility class)
    WindowFunctions() = delete;
 
    static void rectangular(T* output, int size) noexcept
    {
        assert(output != nullptr);
        if (size <= 0) return;
        std::fill(output, output + size, T(1));
    }
 
    static void triangular(T* output, int size, bool periodic = true) noexcept
    {
        assert(output != nullptr);
        if (size <= 0) return;
        if (size == 1) { output[0] = T(1); return; }
 
        const int N = periodic ? size : size - 1;
        const double halfN = static_cast<double>(N) / 2.0;
 
        for (int i = 0; i < size; ++i)
            output[i] = static_cast<T>(1.0 - std::abs((static_cast<double>(i) - halfN) / halfN));
    }
 
    static void hann(T* output, int size, bool periodic = true) noexcept
    {
        static constexpr double kCoeffs[] = { 0.5, -0.5 };
        cosineSum(output, size, periodic, kCoeffs);
    }
 
    static void hamming(T* output, int size, bool periodic = true) noexcept
    {
        static constexpr double kCoeffs[] = { 0.54, -0.46 };
        cosineSum(output, size, periodic, kCoeffs);
    }
 
    static void blackman(T* output, int size, bool periodic = true) noexcept
    {
        static constexpr double kCoeffs[] = { 0.42, -0.5, 0.08 };
        cosineSum(output, size, periodic, kCoeffs);
    }
 
    static void blackmanHarris(T* output, int size, bool periodic = true) noexcept
    {
        // Harris 1978, 4-term minimum side-lobe coefficients.
        static constexpr double kCoeffs[] = { 0.35875, -0.48829, 0.14128, -0.01168 };
        cosineSum(output, size, periodic, kCoeffs);
    }
 
    static void flatTop(T* output, int size, bool periodic = true) noexcept
    {
        // ISO 18431-2 / MATLAB flattopwin coefficients.
        static constexpr double kCoeffs[] = { 0.21557895, -0.41663158, 0.277263158,
                                              -0.083578947, 0.006947368 };
        cosineSum(output, size, periodic, kCoeffs);
    }
 
    static void kaiser(T* output, int size, T beta, bool periodic = true) noexcept
    {
        assert(output != nullptr);
        if (size <= 0) return;
        if (size == 1) { output[0] = T(1); return; }
 
        const int N = periodic ? size : size - 1;
        const double b = static_cast<double>(beta);
        const double invDenominator = 1.0 / besselI0(b);
 
        for (int i = 0; i < size; ++i)
        {
            const double x = 2.0 * static_cast<double>(i) / static_cast<double>(N) - 1.0;
            // Protected against rounding making (1 - x*x) strictly negative.
            const double arg = b * std::sqrt(std::max(0.0, 1.0 - x * x));
            output[i] = static_cast<T>(besselI0(arg) * invDenominator);
        }
    }
 
    // -- Utility methods -------------------------------------------------------
 
    static void apply(T* signal, const T* window, int size) noexcept
    {
        assert(signal != nullptr && window != nullptr);
 
        // No alignment assumption: std::assume_aligned<32> on a pointer that is
        // not actually 32-byte aligned is undefined behaviour, and these buffers
        // come from arbitrary callers (sub-views, vectors -> typically 16-byte).
        for (int i = 0; i < size; ++i)
            signal[i] *= window[i];
    }
 
    [[nodiscard]] static T coherentGain(const T* window, int size) noexcept
    {
        assert(window != nullptr);
        if (size <= 0) return T(0);
 
        double sum = 0.0;
        for (int i = 0; i < size; ++i)
            sum += static_cast<double>(window[i]);
 
        return static_cast<T>(sum / static_cast<double>(size));
    }
 
    [[nodiscard]] static T energyGain(const T* window, int size) noexcept
    {
        assert(window != nullptr);
        if (size <= 0) return T(0);
 
        double sumSq = 0.0;
        for (int i = 0; i < size; ++i)
        {
            const double w = static_cast<double>(window[i]);
            sumSq += w * w;
        }
 
        return static_cast<T>(std::sqrt(sumSq / static_cast<double>(size)));
    }
 
private:
    template <std::size_t NumTerms>
    static void cosineSum(T* output, int size, bool periodic,
                          const double (&coeffs)[NumTerms]) noexcept
    {
        assert(output != nullptr);
        if (size <= 0) return;
        if (size == 1) { output[0] = T(1); return; }
 
        const int N = periodic ? size : size - 1;
        const double twoPiInvN = 2.0 * std::numbers::pi_v<double> / static_cast<double>(N);
 
        for (int i = 0; i < size; ++i)
        {
            const double x = twoPiInvN * static_cast<double>(i);
            double acc = coeffs[0];
            for (std::size_t t = 1; t < NumTerms; ++t)
                acc += coeffs[t] * std::cos(static_cast<double>(t) * x);
            output[i] = static_cast<T>(acc);
        }
    }
 
    [[nodiscard]] static double besselI0(double x) noexcept
    {
        double sum = 1.0;
        double term = 1.0;
        const double halfX = x / 2.0;
 
        for (int k = 1; k < 50; ++k)
        {
            const double f = halfX / static_cast<double>(k);
            term *= f * f;
            sum += term;
            if (term < sum * 1e-12) break; // converged
        }
 
        return sum;
    }
};
 
} // namespace dspark