C++ ADAA tanh() saturator

tanh-saturator.h

// Tanh() saturator with ADAA (1st-order)
// (c) 2026 Signalsmith Audio
// 0BSD Licence, or you can re-distribution with any license you like, just replace this line
#pragma once

#include <cmath>

template<typename Sample>
class TanhSaturator {
	// Accuracy limit to avoid divide-by-(almost)-zero
	static constexpr Sample accuracyLimit = Sample(sizeof(Sample) == sizeof(double) ? 1e-6 : 1e-3);
	Sample prevX = 0, prevIntegralX = 0;

	static Sample curveDirect(Sample x) {
		return std::tanh(x);
	}
	static Sample curveIntegral(Sample x) {
		if (std::abs(x) > 10) return std::abs(x);
		return std::log(2*std::cosh(x));
	}
	
public:
	void reset() {
		prevX = 0;
		prevIntegralX = curveIntegral(prevX);
	}

	inline Sample operator()(Sample x) {
		Sample integralX = curveIntegral(x);
		Sample diffX = x - prevX;
		Sample diffIntegralX = integralX - prevIntegralX;
		prevX = x;
		prevIntegralX = integralX;
		
		if (abs(diffX) > accuracyLimit) {
			return diffIntegralX/diffX;
		} else {
			return curveDirect(x - diffX/2);
		}
	}
};

test.cpp

#include <iostream>
#define LOG_EXPR(expr) std::cout << #expr << " = " << (expr) << std::endl;

#include "./tanh-saturator.h"

#include "signalsmith-util/wav.h" // https://github.com/geraintluff/util

int main() {
	TanhSaturator<float> saturator;

	Wav wav;
	wav.channels = 1;
	wav.sampleRate = 44100;
	float sweepPhase = 0;
	float oscPhase = 0;
	
	while (sweepPhase < 8) {
		// Increase the sine gain exponentially, driving harder into the clipper
		float gain = std::exp(std::floor(sweepPhase))/10;
		
		float sweepPhaseMod = sweepPhase - std::floor(sweepPhase);
		sweepPhase += 1/float(5*wav.sampleRate); // each sweep takes 5s
		float freq = 0.25*sweepPhaseMod*sweepPhaseMod; // Goes up to Nyquist/2, to mimic 2x oversampling
		
		oscPhase += freq;
		oscPhase -= std::floor(oscPhase);
		
		float x = std::sin(oscPhase*float(2*M_PI))*gain;
		
		wav.samples.push_back(saturator(x));
	}
	
	wav.write("sweep.wav");
}