Dither.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 <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <limits>
#include <type_traits>
 
namespace dspark {
 
template <typename T>
class Dither
{
    static_assert(std::is_floating_point_v<T>, "Dither requires a floating-point sample type");
 
public:
    explicit Dither(int targetBits = 16, bool noiseShaping = false) noexcept
        : noiseShaping_(noiseShaping)
    {
        static std::atomic<uint32_t> seedGenerator{ 0x193A6B54u };
        // Increment by golden ratio to ensure diverse initial states across instances
        rngState_ = seedGenerator.fetch_add(0x9E3779B9u, std::memory_order_relaxed);
        // xorshift32 has a fixed point at 0: if the rolling seed ever lands
        // there, the dither would fall silent for that instance.
        if (rngState_ == 0) rngState_ = 1;
 
        setTargetBitDepth(targetBits);
    }
 
    void setTargetBitDepth(int bits) noexcept
    {
        // float (32-bit IEEE 754) only has 24 bits of mantissa precision.
        constexpr int maxBits = (sizeof(T) == 4) ? 24 : 32;
        targetBits_ = std::clamp(bits, 8, maxBits);
 
        // Levels per unit: 2^(N-1). The representable grid is asymmetric,
        // [-levels, levels - 1], mirroring the signed-integer formats.
        int64_t levels = int64_t(1) << (targetBits_ - 1);
 
        quantScale_ = static_cast<T>(levels);
        quantStep_  = T(1) / quantScale_;
        loLevel_    = -quantScale_;
        hiLevel_    = quantScale_ - T(1);
        // Feedback cap: a legitimate shaping error never exceeds 1.5 LSB
        // (TPDF within (-1, 1] LSB + quantiser within +-0.5 LSB); see
        // processSampleInternal for why the cap is needed at all.
        errorCap_   = T(2) * quantStep_;
    }
 
    void setNoiseShaping(bool enabled) noexcept { noiseShaping_ = enabled; }
 
    void reset() noexcept
    {
        std::fill(errorState_.begin(), errorState_.end(), T(0));
    }
 
    [[nodiscard]] T processSample(T input, int channel = 0) noexcept
    {
        if (noiseShaping_ && channel >= 0 && channel < kMaxChannels)
            return processSampleInternal<true>(input, channel);
        else
            return processSampleInternal<false>(input, channel);
    }
 
    void processBlock(T* data, int numSamples, int channel = 0) noexcept
    {
        if (noiseShaping_ && channel >= 0 && channel < kMaxChannels)
        {
            for (int i = 0; i < numSamples; ++i)
                data[i] = processSampleInternal<true>(data[i], channel);
        }
        else
        {
            for (int i = 0; i < numSamples; ++i)
                data[i] = processSampleInternal<false>(data[i], channel);
        }
    }
 
    [[nodiscard]] int getTargetBitDepth() const noexcept { return targetBits_; }
    [[nodiscard]] T getQuantisationStep() const noexcept { return quantStep_; }
 
private:
    static constexpr int kMaxChannels = 16;
 
    template <bool ApplyNoiseShaping>
    [[nodiscard]] inline T processSampleInternal(T input, int channel) noexcept
    {
        // 2 LSB peak-to-peak TPDF noise
        const T noise = (nextRandom() + nextRandom()) * quantStep_;
 
        T preQuantise = input;
 
        if constexpr (ApplyNoiseShaping)
        {
            preQuantise -= errorState_[static_cast<size_t>(channel)];
        }
 
        // Add dither to the signal (w[n] = v[n] + d[n])
        const T dithered = preQuantise + noise;
 
        // Quantise on the integer level grid and clamp to the representable
        // range of the target depth. Ties are broken by the dither, so the
        // rounding mode of nearbyint (default round-to-nearest) is fine and
        // compiles to a single instruction on SSE4.1/NEON.
        const T level = std::clamp(std::nearbyint(dithered * quantScale_), loLevel_, hiLevel_);
        const T quantised = level * quantStep_;
 
        if constexpr (ApplyNoiseShaping)
        {
            // Feed back the TOTAL requantisation error e[n] = y[n] - v[n]
            // (dither + quantiser), the classic error-feedback structure
            // with the dither inside the loop: the whole noise floor is
            // shaped by (1 - z^-1), not just the quantiser error.
            // The cap bounds the loop when the level clamp engages (input at
            // or beyond positive full scale): without it the clip error
            // accumulates sample after sample and pins the output at full
            // scale long after the overload has passed.
            const T err = quantised - preQuantise;
            errorState_[static_cast<size_t>(channel)] = std::clamp(err, -errorCap_, errorCap_);
        }
 
        return quantised;
    }
 
    [[nodiscard]] inline T nextRandom() noexcept
    {
        rngState_ ^= rngState_ << 13;
        rngState_ ^= rngState_ >> 17;
        rngState_ ^= rngState_ << 5;
 
        constexpr T scale = T(1) / static_cast<T>(std::numeric_limits<uint32_t>::max());
        return static_cast<T>(rngState_) * scale - T(0.5);
    }
 
    int targetBits_ = 16;
    T quantStep_ = T(1) / T(32768);
    T quantScale_ = T(32768);
    T loLevel_ = T(-32768);
    T hiLevel_ = T(32767);
    T errorCap_ = T(2) / T(32768);
    bool noiseShaping_ = false;
    uint32_t rngState_;
    std::array<T, kMaxChannels> errorState_{};
};
 
} // namespace dspark