Reverb.cpp

// Copyright (C) 2025, 2026 Theo Niessink <theo@taletn.com>
// This work is free. You can redistribute it and/or modify it under the
// terms of the Do What The Fuck You Want To Public License, Version 2,
// as published by Sam Hocevar. See http://www.wtfpl.net/ for more details.

#include "Reverb.h"

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#ifndef NDEBUG
	#include <stdio.h>
	#include "WDL/heapbuf.h"
#endif

namespace TG101 {

// (9.4 MHz)/224
const double Reverb::mNativeSampleRate = 41964.285714285714;

// Delay 1
const unsigned short Reverb::mMaxDelayTbl[9] =
{
	1323,  // mERDelayTbl[6][0]
	12701, // mERDelayTbl[6][2]

	801,   // mDiffDelayTbl[0][0]
	547,   // mDiffDelayTbl[0][1]
	722,   // mDiffDelayTbl[0][2]
	715,   // mDiffDelayTbl[0][3]

	4126,  // mCombDelayTbl[0][0][0]
	3728,  // mCombDelayTbl[0][1][0]
	3192   // mCombDelayTbl[0][2][0]
};

const unsigned short Reverb::mERDelayTbl[8][3] =
{
	// ER delay   Tap 1  Tap 2
	{ 4521,       3969,  4455  }, // Hall 1
	{ 7167,       6615,  7101  }, // Hall 2
	{ 1875,       1323,  1809  }, // Room 1
	{ 993,        441,   927   }, // Room 2
	{ 666,        396,   617   }, // Plate 1
	{ 2916,       2646,  2867  }, // Plate 2

	// Pre-delay  Tap L  Tap R
	{ 1323,       9570,  12701 }, // Delay 1
	{ 1323,       6703,  2601  }  // Delay 2
};

const float Reverb::mERGainTbl[3][4] =
{
	// Tap 1 (L)   Tap 1 (R)   Tap 2 (L)    Tap 2 (R)
	{ 0.0625f,     0.0625f,    0.078125f,   0.08203125f }, // Hall
	{ 0.046875f,   0.046875f,  0.0625f,     0.06640625f }, // Room
	{ 0.0859375f,  0.0859375f, 0.08984375f, 0.09375f    }  // Plate
};

const float Reverb::mDelayFeedbackTbl[2][2] =
{
	// Feedback (L)  Feedback (R)
	{ 0.3671875f,    0.375f     }, // Delay 1
	{ 0.3359375f,    0.3359375f }  // Delay 2
};

/* static double LowPassCutoff(const int idx)
{
	assert(idx >= 0 && idx < 2);

	const double g = tan(M_PI * 8000.0 / 44100.0);
	const int cutoff = (int)(g / (1.0 + g) * 65536.0 + 0.5);

	const int coeff = (32768 - cutoff + 4) >> (3 + 5);
	return atan((double)(128 - coeff) / (double)(128 - idx + coeff));
} */

const double Reverb::mLowPassCutoffTbl[2] =
{
	0.57004046369971026, // Hall/Room/Plate
	0.57296614288870616  // Delay
};

const unsigned short Reverb::mDiffDelayTbl[3][4] =
{
	// Input     Output
	{ 801, 547,  722, 715 }, // Hall
	{ 772, 539,  722, 715 }, // Room
	{ 629, 676,  371, 386 }  // Plate
};

const float Reverb::mDiffCoeffTbl[4][4] =
{
	// Feedback   Feedforward    Feedback    Feedforward
	{ 0.3046875f, -0.703125f,    0.6875f,    1.6875f     }, // Hall/Room input
	{ 0.2734375f, -0.734375f,    0.71875f,   1.71875f    }, // Hall/Room input
	{ 0.296875f,  -0.703125f,    0.6953125f, 1.6875f     }, // Plate input
	{ 0.59375f,   -0.439453125f, 0.6875f,    0.52734375f }  // Hall/Room/Plate output
};

/* static double CombDamp(const int type)
{
	assert(type >= 0 && type <= 5);

	const int coeff = (32768 - mDampCutoffTbl[type] + 4) >> (3 + 5);
	return atan((double)(128 - coeff) / (double)(128 + coeff));
} */

const double Reverb::mCombDampTbl[6] =
{
	0.60760793440477224, // Hall 1
	0.71520119232557779, // Hall 2
	0.57004046369971026, // Room 1
	0.43349433798248352, // Room 2
	0.57004046369971026, // Plate 1
	0.64568928910828466  // Plate 2
};

const unsigned short Reverb::mCombDelayTbl[3][4][3] =
{
	// Hall
	{
		// Delay  Tap L  Tap R
		{ 4126,   2451,  2387 },
		{ 3728,   1735,  1799 },
		{ 3192,   1243,  1346 },
		{ 2530,   839,   876  }
	},

	// Room
	{
		{ 2647,  2338,  2341  },
		{ 2235,  2015,  1940  },
		{ 1910,  1646,  1740  },
		{ 1658,  1054,  1055  }
	},

	// Plate
	{
		{ 1968,  1157,  1198  },
		{ 1658,  986,   921   },
		{ 1207,  697,   740   },
		{ 1079,  391,   592   }
	}
};

/* static int DampCutoff(const int type)
{
	assert(type >= 0 && type <= 5);

	static const short freq[6] = { 8500, 10000, 8000, 6000, 8000, 9000 };
	const double g = tan(M_PI * (double)freq[type] / 44100.0);

	const double coeff = g / (1.0 + g);
	return (int)(coeff * 65536.0 + 0.5);
} */

const unsigned short Reverb::mDampCutoffTbl[6] =
{
	26809, // Hall 1
	30371, // Hall 2
	25595, // Room 1
	20510, // Room 2
	25595, // Plate 1
	28008  // Plate 2
};

/* static double CombFeedback(const int type)
{
	assert(type >= 0 && type <= 5);

	const int coeff = (32768 - mDampCutoffTbl[type] + 4) >> (3 + 5);
	return 2.0 / (double)(128 - coeff);
} */

const float Reverb::mCombFeedbackTbl[6] =
{
	0.0190476190f, // Hall 1
	0.0168067227f, // Hall 2
	0.02f,         // Room 1
	0.0246913580f, // Room 2
	0.02f,         // Plate 1
	0.0181818182f  // Plate 2
};

const float Reverb::mDelayMixTbl[2][4] =
{
	// Delay (L)     Delay (R)      Reverb (L)  Reverb (R)
	{ 0.143402100f,  0.143402100f,  0.33125f,   0.33125f }, // Delay 1
	{ 0.0484466553f, 0.0484466553f, 0.6375f,    0.6375f  }  // Delay 2
};

const unsigned char Reverb::mOutputLevelTbl[13] =
{
	10,  // -6 dB
	12,  // -5
	13,  // -4
	14,  // -3
	16,  // -2
	18,  // -1

	64,  //  0 dB
	71,  // +1
	80,  // +2
	90,  // +3
	101, // +4
	113, // +5
	127  // +6
};

static int DelayLength(int len, const double ratio)
{
	len = (int)((double)len * ratio + 0.5);
	return wdl_max(len, 1);
}

Reverb::Reverb()
{
	mReverbType = -1;

	mBufSize = 0;
	mBuf = NULL;

	mOutputLevel = 0.0f;
}

Reverb::~Reverb()
{
	free(mBuf);
}

void Reverb::SetSampleRate(const double sampleRate)
{
	mReverbType = -1;

	const double ratio = sampleRate / mNativeSampleRate;
	mSampleRatio = ratio;

	int size = DelayLength(4, ratio);
	mLatency.mLen = size;

	for (int i = 0; i < 2; ++i)
	{
		mLowPassCoeff[i] = GetLowPassCoeff(mLowPassCutoffTbl[i], ratio);
	}

	for (int i = 0; i < 9; ++i)
	{
		size += DelayLength(mMaxDelayTbl[i], ratio);
	}

	free(mBuf);
	mBuf = (float*)malloc((unsigned int)size * sizeof(float));
}

void Reverb::SetReverbType(const EReverbType type)
{
	assert(type >= 0 && type < kNumReverbTypes);

	if (!mBuf) return;
	mReverbType = type;

	mLatency.mIdx = 0;
	int size = mLatency.mLen;

	const double ratio = mSampleRatio;

	const unsigned short* const delay = mERDelayTbl[(unsigned int)type];
	int len[3];

	for (int i = 0; i < 3; ++i)
	{
		len[i] = DelayLength(delay[i], ratio);
	}

	mDelay.mIdx = 0;
	size += len[0];

	if (type < kReverbDelay1)
	{
		mDelay.mLen = len[0];

		mDelayTap[0] = len[1];
		mDelayTap[1] = len[2];
	}
	else
	{
		size += mDelay.mLen = wdl_max(len[1], len[2]);

		mDelayTap[0] = len[1];
		mDelayTap[1] = len[2];

		mPreDelay.mIdx = 0;
		mPreDelay.mLen = len[0];
	}

	mLowPassState = 0.0;

	size += SetDiffuserType(type, ratio);
	size += SetCombFilterType(type, ratio);

	mBufSize = size;
	memset(mBuf, 0, (unsigned int)size * sizeof(float));
}

int Reverb::SetDiffuserType(int type, const double ratio)
{
	assert(type >= 0 /* && type < kNumReverbTypes */);

	type = type < kReverbDelay1 ? type : kReverbHall1;
	const unsigned short* delay = mDiffDelayTbl[(unsigned int)type >> 1];

	DelayLine* diff = mInputDiff;
	int size = 0;

	for (int i = 0; i < 2; ++i)
	{
		for (int j = 0; j < 2; ++j)
		{
			diff[j].mIdx = 0;
			size += diff[j].mLen = DelayLength(delay[j], ratio);
		}

		delay += 2;
		diff = mOutputDiff;
	}

	return size;
}

int Reverb::SetCombFilterType(int type, const double ratio)
{
	assert(type >= 0 /* && type < kNumReverbTypes */);

	const int n = type < kReverbDelay1 ? 4 : 3;
	type = type < kReverbDelay1 ? type : kReverbHall1;

	mHighDamp = GetLowPassCoeff(mCombDampTbl[(unsigned int)type], ratio);
	int size = 0;

	for (int i = 0; i < n; ++i)
	{
		const unsigned short* const delay = mCombDelayTbl[(unsigned int)type >> 1][i];
		int len[3];

		for (int j = 0; j < 3; ++j)
		{
			len[j] = DelayLength(delay[j], ratio);
		}

		mComb[i].mDelay.mIdx = 0;
		size += mComb[i].mDelay.mLen = len[0];

		mComb[i].mTap[0] = len[1];
		mComb[i].mTap[1] = len[2];

		#ifndef NDEBUG
		mComb[i].mGain = mComb[i].mFeedback = 0.0f;
		#endif

		mComb[i].mFilterState = 0.0;
	}

	return size;
}

double Reverb::GetLowPassCoeff(const double cutoff, const double ratio)
{
	double x = cutoff / ratio;
	x = wdl_min(x, 1.5550883635269477); // 0.495*pi

	// g = tan(x)
	const double y1 = sin(x);
	const double y2 = cos(x);

	// g/(1 + g)
	return y1 / (y1 + y2);
}

void Reverb::SetReverbFeedback(const unsigned short coeff[4][2], const int num)
{
	assert(num == 3 || num == 4);

	const int type = mReverbType;
	if (type < 0) return;

	const unsigned int reverb = type < kReverbDelay1 ? type : kReverbHall1;

	const int mul = mDampCutoffTbl[reverb];
	const float div = mCombFeedbackTbl[reverb];

	static const int scale4 = 32768; // 65536*sqrt(1/4)
	static const int scale3 = 37837; // 65536*sqrt(1/3)

	const int scale = num == 4 ? scale4 : scale3;

	for (int i = 0; i < num; ++i)
	{
		int feedback = coeff[i][0];
		int gain = coeff[i][1];

		assert(feedback >= 3 && feedback <= 4095);
		assert(gain >= 89 && gain <= 4096);

		feedback = (unsigned int)(feedback * mul) >> 21;

		gain = (unsigned int)(gain * scale) >> 21;
		gain = wdl_max(gain, 3);

		mComb[i].mFeedback = (float)feedback * div;
		mComb[i].mGain = (float)gain * 0.0078125f; // 1/128
	}

	if (type >= kReverbDelay1) mComb[0].mFeedback *= 0.9921875f;
}

void Reverb::SetOutputLevel(int level)
{
	assert(level >= 24 && level <= 70);

	unsigned int idx;
	int scale;

	if (level < 64)
	{
		level -= 24;

		scale = (unsigned int)level / 6;
		idx = (unsigned int)level % 6;
	}
	else
	{
		level -= 64;

		scale = 4;
		idx = level + 6;
	}

	level = mOutputLevelTbl[idx] << scale;
	mOutputLevel = (float)level * 0.0009765625f; // 1/1024
}

void Reverb::Reset()
{
	const int type = mReverbType;
	const unsigned int size = mBufSize;

	mLatency.mIdx = 0;
	mDelay.mIdx = 0;
	mPreDelay.mIdx = 0;

	mLowPassState = 0.0;

	mInputDiff[0].mIdx = 0;
	mInputDiff[1].mIdx = 0;

	for (int i = 0; i < 4; ++i)
	{
		mComb[i].mDelay.mIdx = 0;
		mComb[i].mFilterState = 0.0;
	}

	mOutputDiff[0].mIdx = 0;
	mOutputDiff[1].mIdx = 0;

	if (type >= 0) memset(mBuf, 0, size * sizeof(float));
}

void Reverb::Process(float input, float output[2])
{
	const int type = mReverbType;
	float* buf = mBuf;

	if (type >= 0)
	{
		buf = ProcessDelayLine(input, &input, &mLatency, buf);

		float reverb[2];

		if (type < kReverbDelay1)
			buf = ProcessReverb(input, reverb, type, buf);
		else
			buf = ProcessDelay(input, reverb, type, buf);

		assert(buf == mBuf + mBufSize);

		const float level = mOutputLevel;

		reverb[0] *= level;
		reverb[1] *= level;

		output[0] = reverb[0];
		output[1] = reverb[1];
	}
	else
	{
		memset(output, 0, 2 * sizeof(float));
	}
}

float* Reverb::ProcessDelayLine(const float input, float* const output, DelayLine* const delay, float* const buf)
{
	unsigned int idx = delay->mIdx;
	const unsigned int len = delay->mLen;

	*output = buf[idx];
	buf[idx] = input;

	delay->mIdx = ++idx < len ? idx : 0;

	return buf + len;
}

float* Reverb::ProcessReverb(float input, float output[2], const int type, float* buf)
{
	assert(type >= kReverbHall1 && type <= kReverbPlate2);

	float delay[3];
	buf = ProcessEarlyReflect(input, delay, type, buf);

	static const float gain = 0.98f;
	const float preDelay = delay[2];

	input = ProcessLowPassFilter(preDelay * gain, &mLowPassState, mLowPassCoeff[0]);

	const unsigned int plate = type & kReverbPlate1;

	const float* inCoeff = mDiffCoeffTbl[plate >> 1];
	const unsigned int delta = !plate ? 4 : 0;

	for (int i = 0; i < 2; ++i)
	{
		buf = ProcessDiffuser(input, &input, &mInputDiff[i], buf, inCoeff);
		inCoeff += delta;
	}

	const double damp = mHighDamp;
	float reverb[2] = { 0.0f, 0.0f };

	for (int i = 0; i < 4; ++i)
	{
		buf = ProcessCombFilter(input, reverb, &mComb[i], buf, damp);
	}

	const float* const outCoeff = mDiffCoeffTbl[3];

	for (int i = 0; i < 2; ++i)
	{
		buf = ProcessDiffuser(reverb[i], &reverb[i], &mOutputDiff[i], buf, outCoeff);
	}

	reverb[0] += delay[0];
	reverb[1] += delay[1];

	output[0] = reverb[0];
	output[1] = reverb[1];

	return buf;
}

float* Reverb::ProcessDelay(float input, float output[2], const int type, float* buf)
{
	assert(type == kReverbDelay1 || type == kReverbDelay2);

	const float* const feedback = mDelayFeedbackTbl[type & 1];

	float delay[4];
	buf = ProcessDelayFeedback(input, delay, feedback, buf);

	float preDelay = input;

	if (type == kReverbDelay2)
	{
		static const float dry = 0.25f;
		static const float wet = 0.7421875f;

		input *= dry;

		delay[0] *= wet;
		delay[1] *= wet;

		delay[0] += input;
		delay[1] += input;

		preDelay = (delay[0] + delay[1]) * 0.5f;
	}

	buf = ProcessDelayLine(preDelay, &preDelay, &mPreDelay, buf);

	static const float gain = 0.974018555f;
	input = ProcessLowPassFilter(preDelay * gain, &mLowPassState, mLowPassCoeff[1]);

	const float* inCoeff = mDiffCoeffTbl[0]; // kReverbHall1
	static const unsigned int delta = 4;

	for (int i = 0; i < 2; ++i)
	{
		buf = ProcessDiffuser(input, &input, &mInputDiff[i], buf, inCoeff);
		inCoeff += delta;
	}

	float* const reverb = &delay[2];

	const double damp = mHighDamp;
	reverb[1] = reverb[0] = 0.0f;

	for (int i = 0; i < 3; ++i)
	{
		buf = ProcessCombFilter(input, reverb, &mComb[i], buf, damp);
	}

	const float* const outCoeff = mDiffCoeffTbl[3];

	for (int i = 0; i < 2; ++i)
	{
		buf = ProcessDiffuser(reverb[i], &reverb[i], &mOutputDiff[i], buf, outCoeff);
	}

	const unsigned int mix = type & 1;

	for (int i = 0; i < 4; ++i)
	{
		delay[i] *= mDelayMixTbl[mix][i];
	}

	delay[0] += delay[2];
	delay[1] += delay[3];

	output[0] = delay[0];
	output[1] = delay[1];

	return buf;
}

float* Reverb::ProcessEarlyReflect(const float input, float output[3], const int type, float* const buf)
{
	assert(type >= kReverbHall1 && type <= kReverbPlate2);

	unsigned int idx = mDelay.mIdx;
	const unsigned int len = mDelay.mLen;

	unsigned int ofs[2] = { idx, idx };

	ofs[0] -= mDelayTap[0];
	ofs[1] -= mDelayTap[1];

	ofs[0] += (int)ofs[0] < 0 ? len : 0;
	ofs[1] += (int)ofs[1] < 0 ? len : 0;

	const float preDelay = buf[ofs[0]];
	float tap[4] = { preDelay, preDelay, buf[ofs[1]], buf[idx] };

	buf[idx] = input;
	mDelay.mIdx = ++idx < len ? idx : 0;

	const unsigned int er = (unsigned int)type >> 1;

	for (int i = 0; i < 4; ++i)
	{
		tap[i] *= mERGainTbl[er][i];
	}

	tap[0] += tap[2];
	tap[1] += tap[3];

	output[0] = tap[0];
	output[1] = tap[1];
	output[2] = preDelay;

	return buf + len;
}

float* Reverb::ProcessDelayFeedback(const float input, float output[2], const float coeff[2], float* const buf)
{
	unsigned int idx = mDelay.mIdx;
	const unsigned int len = mDelay.mLen;

	unsigned int ofs[2] = { idx, idx };

	ofs[0] -= mDelayTap[0];
	ofs[1] -= mDelayTap[1];

	ofs[0] += (int)ofs[0] < 0 ? len : 0;
	ofs[1] += (int)ofs[1] < 0 ? len : 0;

	float tap[2];

	tap[0] = buf[ofs[0]];
	tap[1] = buf[ofs[1]];

	float feedback[2] = { coeff[0], coeff[1] };

	feedback[0] *= tap[0];
	feedback[1] *= tap[1];

	feedback[0] += feedback[1];
	buf[idx] = input + feedback[0];

	mDelay.mIdx = ++idx < len ? idx : 0;

	output[0] = tap[0];
	output[1] = tap[1];

	return buf + len;
}

float Reverb::ProcessLowPassFilter(const float input, double* const state, const double coeff)
{
	const double x = (double)input;
	const double s = *state;

	const double v = (x - s) * coeff;
	const double y = v + s;
	*state = y + v;

	return (float)y;
}

float* Reverb::ProcessDiffuser(const float input, float* const output, DelayLine* const diff, float* const buf, const float coeff[4])
{
	unsigned int idx = diff->mIdx;
	const unsigned int len = diff->mLen;

	const float& x = input;
	const float y = buf[idx];

	const float feedback = x * coeff[0] + y * coeff[2];
	const float feedfwd = x * coeff[1] + y * coeff[3];

	buf[idx] = feedback;
	*output = feedfwd;

	diff->mIdx = ++idx < len ? idx : 0;

	return buf + len;
}

float* Reverb::ProcessCombFilter(const float input, float output[2], CombFilter* const comb, float* const buf, const double damp)
{
	unsigned int idx = comb->mDelay.mIdx;
	const unsigned int len = comb->mDelay.mLen;

	unsigned int ofs[2] = { idx, idx };

	ofs[0] -= comb->mTap[0];
	ofs[1] -= comb->mTap[1];

	ofs[0] += (int)ofs[0] < 0 ? len : 0;
	ofs[1] += (int)ofs[1] < 0 ? len : 0;

	output[0] += buf[ofs[0]];
	output[1] += buf[ofs[1]];

	const float feedback = comb->mFeedback;
	const float gain = comb->mGain;

	const float x = buf[idx] * feedback;
	const float y = ProcessLowPassFilter(x, &comb->mFilterState, damp);

	buf[idx] = input * gain + y;
	comb->mDelay.mIdx = ++idx < len ? idx : 0;

	return buf + len;
}

#ifndef NDEBUG

void Reverb::ExtractIR(const char* const filename, int maxLen)
{
	assert(mReverbType >= 0);

	if (maxLen <= 0)
	{
		const double sampleRate = mNativeSampleRate * mSampleRatio;
		maxLen = (int)ceil(20.0 * sampleRate);
	}

	WDL_TypedBuf<float> buf;
	float* ptr = buf.ResizeOK(maxLen * 2);

	if (!ptr) return;

	Reset();
	const float level = mOutputLevel;

	mOutputLevel = 1.0f;
	float input = 1.0f;

	for (int i = 0; i < maxLen; ++i)
	{
		Process(input, ptr);
		input = 0.0f;
		ptr += 2;
	}

	Reset();
	mOutputLevel = level;

	// "Silence will fall"
	bool silence = true;

	float* end = ptr;
	ptr = buf.GetFast();

	static const float threshold = 7.62939453e-06f; // 2^(-17)

	while (silence && end > ptr)
	{
		end -= 2;
		silence = fabs(end[0]) < threshold && fabs(end[1]) < threshold;
	}

	FILE* const fp = fopen(filename, "wb");

	if (fp)
	{
		fwrite(ptr, sizeof(float), end + 2 - ptr, fp);
		fclose(fp);
	}
}

#endif // NDEBUG

} // namespace TG101