TruePeakDetector.h

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include "DspMath.h"
#include "SimdOps.h"
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
 
namespace dspark {
 
template <FloatType T, int MaxChannels = 16>
class TruePeakDetector
{
public:
    void reset() noexcept
    {
        for (auto& s : states_)
            s = {};
    }
 
    [[nodiscard]] T processSample(T sample, int channel) noexcept
    {
        assert(channel >= 0 && channel < MaxChannels);
        // Release-safe clamp to the nearest valid channel. (Redirecting an
        // out-of-range index to channel 0 would corrupt that channel's
        // history with a foreign stream.)
        channel = std::clamp(channel, 0, MaxChannels - 1);
        return processState(sample, states_[static_cast<size_t>(channel)]);
    }
 
    [[nodiscard]] T getTailPeak(int channel) const noexcept
    {
        assert(channel >= 0 && channel < MaxChannels);
        channel = std::clamp(channel, 0, MaxChannels - 1);
        auto state = states_[static_cast<size_t>(channel)];
        T peak = T(0);
        for (int i = 0; i < kTaps - 1; ++i)
            peak = std::max(peak, processState(T(0), state));
        return peak;
    }
 
    [[nodiscard]] static constexpr int getLatency() noexcept { return kTaps / 2; }
 
    [[nodiscard]] static constexpr int getTaps() noexcept { return kTaps; }
 
private:
    static constexpr int kTaps     = 12;  
    static constexpr int kPhases   = 4;   
    static constexpr int kHistSize = 16;  
    static constexpr int kHistMask = kHistSize - 1;
 
    struct State
    {
        T history[kHistSize * 2] = {};
        int writePos = 0;
    };
 
    [[nodiscard]] static T processState(T sample, State& tp) noexcept
    {
        // Mirroring makes the latest FIR window contiguous for the shared dot product.
        tp.history[static_cast<size_t>(tp.writePos)] = sample;
        tp.history[static_cast<size_t>(tp.writePos + kHistSize)] = sample;
        const int newest = tp.writePos;
        tp.writePos = (tp.writePos + 1) & kHistMask;
        const T* window = &tp.history[static_cast<size_t>(newest + kHistSize - (kTaps - 1))];
        T peak = std::abs(sample);
        for (int phase = 0; phase < kPhases; ++phase)
        {
            const T interp = simd::dotProduct(
                kReversedCoeffs[static_cast<size_t>(phase)].data(), window, kTaps);
            const T a = std::abs(interp);
            if (a > peak) peak = a;
        }
        return peak;
    }
 
    [[nodiscard]] static constexpr std::array<std::array<T, kTaps>, kPhases>
    buildReversedPhaseTable() noexcept
    {
        constexpr double kAnnex2[kPhases][kTaps] = {
            {  0.0017089843750,  0.0109863281250, -0.0196533203125,
               0.0332031250000, -0.0594482421875,  0.1373291015625,
               0.9721679687500, -0.1022949218750,  0.0476074218750,
              -0.0266113281250,  0.0148925781250, -0.0083007812500 },
            { -0.0291748046875,  0.0292968750000, -0.0517578125000,
               0.0891113281250, -0.1665039062500,  0.4650878906250,
               0.7797851562500, -0.2003173828125,  0.1015625000000,
              -0.0582275390625,  0.0330810546875, -0.0189208984375 },
            { -0.0189208984375,  0.0330810546875, -0.0582275390625,
               0.1015625000000, -0.2003173828125,  0.7797851562500,
               0.4650878906250, -0.1665039062500,  0.0891113281250,
              -0.0517578125000,  0.0292968750000, -0.0291748046875 },
            { -0.0083007812500,  0.0148925781250, -0.0266113281250,
               0.0476074218750, -0.1022949218750,  0.9721679687500,
               0.1373291015625, -0.0594482421875,  0.0332031250000,
              -0.0196533203125,  0.0109863281250,  0.0017089843750 },
        };
 
        std::array<std::array<T, kTaps>, kPhases> result {};
        for (int phase = 0; phase < kPhases; ++phase)
            for (int k = 0; k < kTaps; ++k)
                result[static_cast<size_t>(phase)][static_cast<size_t>(kTaps - 1 - k)]
                    = static_cast<T>(kAnnex2[phase][k]);
        return result;
    }
 
    static constexpr std::array<std::array<T, kTaps>, kPhases> kReversedCoeffs =
        buildReversedPhaseTable();
 
    std::array<State, MaxChannels> states_ {};
};
 
} // namespace dspark