reverb.c

//////////////////
// This is a slightly modified version of Sean Connelly (@voidqk) code that
// uses jni types.
//////////////////

// (c) Copyright 2016, Sean Connelly (@voidqk), http://syntheti.cc
// MIT License
// Project Home: https://github.com/voidqk/sndfilter

#include "reverb.h"
#include <math.h>
#include <stdint.h>

#ifndef M_PI
	#define M_PI 3.14159265358979323846
#endif
#ifndef M_LN2
	#define M_LN2  0.69314718055994530941723212145817656807550013436
#endif

// utility functions
static inline jfloat db2lin(jfloat db){ // dB to linear
	return powf(10.0f, 0.05f * db);
}

static inline jint clampi(jint v, jint min, jint max){
	return v < min ? min : (v > max ? max : v);
}

static inline jfloat clampf(jfloat v, jfloat min, jfloat max){
	return v < min ? min : (v > max ? max : v);
}

static jboolean isprime(jint v){
	if (v < 0)
		return isprime(-v);
	if (v < 2)
		return false;
	if (v <= 3)
		return true;
	if ((v % 2) == 0 || (v % 3) == 0)
		return false;
	jint max = (jint)sqrt((double)v);
	for (jint i = 5; i <= max; i += 6){
		if ((v % i) == 0 || (v % (i + 2)) == 0)
			return false;
	}
	return true;
}

static inline jint nextprime(jint v){
	while (!isprime(v))
		v++;
	return v;
}

// generate a random jfloat [0, 1) using a simple (but good quality) RNG
static inline jfloat randjfloat(){
	static uint32_t seed = 123; // doesn't matter
	static uint32_t i = 456; // doesn't matter
	uint32_t m = 0x5bd1e995;
	uint32_t k = i++ * m;
	seed = (k ^ (k >> 24) ^ (seed * m)) * m;
	uint32_t R = (seed ^ (seed >> 13)) & 0x007FFFFF; // get 23 random bits
	union { uint32_t i; float f; } u = { .i = 0x3F800000 | R };
	return u.f - 1.0;
}

//
//
// component implementation
//
//

// components are in the basic format of `<component>_make` to initialize a structure and
// `<component>_step` to perform a single step with the component

//
// delay
//
static inline void delay_make(sf_rv_delay_st *delay, jint size){
	delay->pos = 0;
	delay->size = clampi(size, 1, SF_REVERB_DS);
	// memset(delay->buf, 0, sizeof(jfloat) * delay->size);
}

static inline jfloat delay_step(sf_rv_delay_st *delay, jfloat v){
	jfloat out = delay->buf[delay->pos];
	delay->buf[delay->pos] = v;
	delay->pos = (delay->pos + 1) % delay->size;
	return out;
}

// delay_get(d, 1) returns the last written value
// delay_get(d, 2) returns the second-last written value
// ..etc
static inline jfloat delay_get(sf_rv_delay_st *delay, jint offset){
	if (offset > delay->size)
		return delay->buf[delay->pos];
	else if (offset <= 0)
		offset = 1;
	jint pos = delay->pos - offset;
	if (pos < 0)
		pos += delay->size;
	return delay->buf[pos];
}

static inline jfloat delay_getlast(sf_rv_delay_st *delay){
	return delay->buf[delay->pos];
}

//
// iir1
//
static inline void iir1_makeLPF(sf_rv_iir1_st *iir1, jint rate, jfloat freq){
	// 1st order IIR lowpass filter (Butterworth)
	freq = clampf(freq, 0, rate / 2);
	jfloat omega2 = (jfloat)M_PI * freq / (jfloat)rate;
	jfloat tano2 = tanf(omega2);
	iir1->b1 = iir1->b2 = tano2 / (1.0f + tano2);
	iir1->a2 = (1.0f - tano2) / (1.0f + tano2);
	iir1->y1 = 0;
}

static inline void iir1_makeHPF(sf_rv_iir1_st *iir1, jint rate, jfloat freq){
	// 1st order IIR highpass filter (Butterworth)
	freq = clampf(freq, 0, rate / 2);
	jfloat omega2 = (jfloat)M_PI * freq / (jfloat)rate;
	jfloat tano2 = tanf(omega2);
	iir1->b1 = 1.0f / (1.0f + tano2);
	iir1->b2 = -iir1->b1;
	iir1->a2 = (1.0f - tano2) / (1.0f + tano2);
	iir1->y1 = 0;
}

static inline jfloat iir1_step(sf_rv_iir1_st *iir1, jfloat v){
	jfloat out = v * iir1->b1 + iir1->y1;
	iir1->y1 = out * iir1->a2 + v * iir1->b2;
	return out;
}

//
// biquad
//
static inline void biquad_makeLPF(sf_rv_biquad_st *biquad, jint rate, jfloat freq, jfloat bw){
	freq = clampf(freq, 0, rate / 2);
	jfloat omega = 2.0f * (jfloat)M_PI * freq / (jfloat)rate;
	jfloat cs = cosf(omega);
	jfloat sn = sinf(omega);
	jfloat alpha = sn * sinhf((jfloat)M_LN2 * 0.5f * bw * omega / sn);
	jfloat a0inv = 1.0f / (1.0f + alpha);
	biquad->b0 = a0inv * (1.0f - cs) * 0.5f;
	biquad->b1 = 2.0f * biquad->b0;
	biquad->b2 = biquad->b0;
	biquad->a1 = a0inv * -2.0f * cs;
	biquad->a2 = a0inv * (1.0f - alpha);
	biquad->xn1 = 0;
	biquad->xn2 = 0;
	biquad->yn1 = 0;
	biquad->yn2 = 0;
}

static inline void biquad_makeLPFQ(sf_rv_biquad_st *biquad, jint rate, jfloat freq, jfloat bw){
	freq = clampf(freq, 0, rate / 2);
	jfloat omega = 2.0f * (jfloat)M_PI * freq / (jfloat)rate;
	jfloat cs = cosf(omega);
	jfloat alpha = sinf(omega) * 2.0f * bw; // different alpha calculation than makeLPF above
	jfloat a0inv = 1.0f / (1.0f + alpha);
	biquad->b0 = a0inv * (1.0f - cs) * 0.5f;
	biquad->b1 = 2.0f * biquad->b0;
	biquad->b2 = biquad->b0;
	biquad->a1 = a0inv * -2.0f * cs;
	biquad->a2 = a0inv * (1.0f - alpha);
	biquad->xn1 = 0;
	biquad->xn2 = 0;
	biquad->yn1 = 0;
	biquad->yn2 = 0;
}

static inline void biquad_makeAPF(sf_rv_biquad_st *biquad, jint rate, jfloat freq, jfloat bw){
	freq = clampf(freq, 0, rate / 2);
	jfloat omega = 2.0f * (jfloat)M_PI * freq / (jfloat)rate;
	jfloat sn = sinf(omega);
	jfloat alpha = sn * sinhf((jfloat)M_LN2 * 0.5f * bw * omega / sn);
	jfloat a0inv = 1.0f / (1.0f + alpha);
	biquad->b0 = a0inv * (1.0f - alpha);
	biquad->b1 = a0inv * -2.0f * cosf(omega);
	biquad->b2 = a0inv * (1.0f + alpha);
	biquad->a1 = biquad->b1;
	biquad->a2 = biquad->b0;
	biquad->xn1 = 0;
	biquad->xn2 = 0;
	biquad->yn1 = 0;
	biquad->yn2 = 0;
}

static inline jfloat biquad_step(sf_rv_biquad_st *biquad, jfloat v){
	jfloat out = v * biquad->b0 + biquad->xn1 * biquad->b1 + biquad->xn2 * biquad->b2 -
		biquad->yn1 * biquad->a1 - biquad->yn2 * biquad->a2;
	biquad->xn2 = biquad->xn1;
	biquad->xn1 = v;
	biquad->yn2 = biquad->yn1;
	biquad->yn1 = out;
	return out;
}

//
// earlyref
//
static inline void earlyref_make(sf_rv_earlyref_st *earlyref, jint rate, jfloat factor, jfloat width){
	static const sf_sample_st delaytbl[18] = {
		// seconds to look backwards
		{ 0.0043f, 0.0053f }, { 0.0215f, 0.0225f }, { 0.0225f, 0.0235f }, { 0.0268f, 0.0278f },
		{ 0.0270f, 0.0290f }, { 0.0298f, 0.0288f }, { 0.0458f, 0.0468f }, { 0.0485f, 0.0475f },
		{ 0.0572f, 0.0582f }, { 0.0587f, 0.0577f }, { 0.0595f, 0.0575f }, { 0.0612f, 0.0622f },
		{ 0.0707f, 0.0697f }, { 0.0708f, 0.0718f }, { 0.0726f, 0.0736f }, { 0.0741f, 0.0751f },
		{ 0.0753f, 0.0763f }, { 0.0797f, 0.0817f }
	};

	earlyref->wet1 = width * 0.5f + 0.5f;
	earlyref->wet2 = (1.0f - width) * 0.5f;

	jint lrdelay = 0.0002f * (jfloat)rate;
	delay_make(&earlyref->delayRL, lrdelay);
	delay_make(&earlyref->delayLR, lrdelay);

	biquad_makeAPF(&earlyref->allpassXL, rate, 740.0f, 4.0f);
	earlyref->allpassXR = earlyref->allpassXL;

	biquad_makeAPF(&earlyref->allpassL, rate, 150.0f, 4.0f);
	earlyref->allpassR = earlyref->allpassL;

	factor *= rate;
	for (jint i = 0; i < 18; i++){
		earlyref->delaytblL[i] = delaytbl[i].L * factor;
		earlyref->delaytblR[i] = delaytbl[i].R * factor;
	}
	delay_make(&earlyref->delayPWL, earlyref->delaytblL[17] + 10);
	delay_make(&earlyref->delayPWR, earlyref->delaytblR[17] + 10);

	iir1_makeLPF(&earlyref->lpfL, rate, 20000.0f);
	earlyref->lpfR = earlyref->lpfL;

	iir1_makeHPF(&earlyref->hpfL, rate, 4.0f);
	earlyref->hpfR = earlyref->hpfL;
}

static inline sf_sample_st earlyref_step(sf_rv_earlyref_st *earlyref, sf_sample_st input){
	static const sf_sample_st gajintbl[18] = {
		{ 0.841f, 0.842f }, { 0.504f, 0.506f }, { 0.491f, 0.489f }, { 0.379f, 0.382f },
		{ 0.380f, 0.300f }, { 0.346f, 0.346f }, { 0.289f, 0.290f }, { 0.272f, 0.271f },
		{ 0.192f, 0.193f }, { 0.193f, 0.192f }, { 0.217f, 0.217f }, { 0.181f, 0.195f },
		{ 0.180f, 0.192f }, { 0.181f, 0.166f }, { 0.176f, 0.186f }, { 0.142f, 0.131f },
		{ 0.167f, 0.168f }, { 0.134f, 0.133f }
	};

	jfloat wetL = 0, wetR = 0;
	delay_step(&earlyref->delayPWL, input.L);
	delay_step(&earlyref->delayPWR, input.R);
	for (jint i = 0; i < 18; i++){
		wetL += gajintbl[i].L * delay_get(&earlyref->delayPWL, earlyref->delaytblL[i]);
		wetR += gajintbl[i].R * delay_get(&earlyref->delayPWR, earlyref->delaytblR[i]);
	}

	jfloat L = delay_step(&earlyref->delayRL, input.R + wetR);
	L = biquad_step(&earlyref->allpassXL, L);
	L = biquad_step(&earlyref->allpassL, earlyref->wet1 * wetL + earlyref->wet2 * L);
	L = iir1_step(&earlyref->hpfL, L);
	L = iir1_step(&earlyref->lpfL, L);

	jfloat R = delay_step(&earlyref->delayLR, input.L + wetL);
	R = biquad_step(&earlyref->allpassXR, R);
	R = biquad_step(&earlyref->allpassR, earlyref->wet1 * wetR + earlyref->wet2 * R);
	R = iir1_step(&earlyref->hpfR, R);
	R = iir1_step(&earlyref->lpfR, R);

	return (sf_sample_st){ L, R };
}

//
// oversample
//
static inline void oversample_make(sf_rv_oversample_st *oversample, jint factor){
	oversample->factor = clampi(factor, 1, SF_REVERB_OF);
	biquad_makeLPFQ(&oversample->lpfU, 2 * oversample->factor, 1.0f,
		0.5773502691896258f); // 1/sqrt(3)
	oversample->lpfD = oversample->lpfU;
}

// output length must be oversample->factor
static inline void oversample_stepup(sf_rv_oversample_st *oversample, jfloat input, jfloat *output){
	if (oversample->factor == 1){
		output[0] = input;
		return;
	}
	output[0] = biquad_step(&oversample->lpfU, input * oversample->factor);
	for (jint i = 1; i < oversample->factor; i++)
		output[i] = biquad_step(&oversample->lpfU, 0);
}

// input length must be oversample->factor
static inline jfloat oversample_stepdown(sf_rv_oversample_st *oversample, jfloat *input){
	if (oversample->factor == 1)
		return input[0];
	for (jint i = 0; i < oversample->factor; i++)
		biquad_step(&oversample->lpfD, input[i]);
	return input[0];
}

//
// dccut
//
static inline void dccut_make(sf_rv_dccut_st *dccut, jint rate, jfloat freq){
	freq = clampf(freq, 0, rate / 2);
	jfloat ang = 2.0f * (jfloat)M_PI * freq / (jfloat)rate;
	jfloat sn = sinf(ang);
	jfloat sqrt3 = 1.7320508075688772f;
	dccut->gain = (sqrt3 - 2.0f * sn) / (sn + sqrt3 * cosf(ang));
	dccut->y1 = 0;
	dccut->y2 = 0;
}

static inline jfloat dccut_step(sf_rv_dccut_st *dccut, jfloat v){
	jfloat out = v - dccut->y1 + dccut->gain * dccut->y2;
	dccut->y1 = v;
	dccut->y2 = out;
	return out;
}

//
// noise
//
static inline void noise_make(sf_rv_noise_st *noise){
	noise->pos = SF_REVERB_NS;
}

static inline jfloat noise_step(sf_rv_noise_st *noise){
	if (noise->pos >= SF_REVERB_NS){
		// need to generate more noise
		noise->pos = 0;
		jint len = SF_REVERB_NS;
		jint tot = 1;
		jfloat r = 0.8f;
		jfloat rmul = 0.7071067811865475f; // 1/sqrt(2)
		noise->buf[0] = 0;
		while (len > 1){
			jfloat left = 0;
			for (jint i = tot - 1; i >= 0; i--){
				jfloat right = left;
				left = noise->buf[i * len];
				jfloat midpojint = (left + right) * 0.5f;
				jfloat newv = midpojint + r * (2.0f * randjfloat() - 1.0f); // displace by random amt
				noise->buf[i * len + (len / 2)] = clampf(newv, -1.0f, 1.0f);
			}
			len /= 2;
			tot *= 2;
			r *= rmul;
		}
	}
	return noise->buf[noise->pos++];
}

//
// lfo
//
static inline void lfo_make(sf_rv_lfo_st *lfo, jint rate, jfloat freq){
	lfo->count = 0;
	lfo->re = 1.0f;
	lfo->im = 0.0f;
	jfloat theta = 2.0f * (jfloat)M_PI * freq / (jfloat)rate;
	lfo->sn = sinf(theta);
	lfo->co = cosf(theta);
}

static inline jfloat lfo_step(sf_rv_lfo_st *lfo){
	jfloat v = lfo->im;
	jfloat re = lfo->re * lfo->co - lfo->im * lfo->sn;
	jfloat im = lfo->re * lfo->sn + lfo->im * lfo->co;
	if (lfo->count++ > 100000){
		// if we've gathered a lot of samples, then it's probably a good idea to make sure our LFO
		// hasn't accumulated a bunch of errors
		lfo->count = 0;
		jfloat leninv = 1.0f / sqrtf(re * re + im * im);
		re *= leninv;
		im *= leninv;
	}
	lfo->re = re;
	lfo->im = im;
	return v;
}

//
// allpass
//
static inline void allpass_make(sf_rv_allpass_st *allpass, jint size, jfloat feedback, jfloat decay){
	allpass->pos = 0;
	allpass->size = clampi(size, 1, SF_REVERB_APS);
	allpass->feedback = feedback;
	allpass->decay = decay;
	// memset(allpass->buf, 0, sizeof(jfloat) * allpass->size);
}

static inline jfloat allpass_step(sf_rv_allpass_st *allpass, jfloat v){
	v += allpass->feedback * allpass->buf[allpass->pos];
	jfloat out = allpass->decay * allpass->buf[allpass->pos] - allpass->feedback * v;
	allpass->buf[allpass->pos] = v;
	allpass->pos = (allpass->pos + 1) % allpass->size;
	return out;
}

//
// allpass2
//
static inline void allpass2_make(sf_rv_allpass2_st *allpass2, jint size1, jint size2, jfloat feedback1,
	jfloat feedback2, jfloat decay1, jfloat decay2){
	allpass2->pos1 = 0;
	allpass2->pos2 = 0;
	allpass2->size1 = clampi(size1, 1, SF_REVERB_AP2S1);
	allpass2->size2 = clampi(size2, 1, SF_REVERB_AP2S2);
	allpass2->feedback1 = feedback1;
	allpass2->feedback2 = feedback2;
	allpass2->decay1 = decay1;
	allpass2->decay2 = decay2;
	// memset(allpass2->buf1, 0, sizeof(jfloat) * allpass2->size1);
	// memset(allpass2->buf2, 0, sizeof(jfloat) * allpass2->size2);
}

static inline jfloat allpass2_step(sf_rv_allpass2_st *allpass2, jfloat v){
	v += allpass2->feedback2 * allpass2->buf2[allpass2->pos2];
	jfloat out = allpass2->decay2 * allpass2->buf2[allpass2->pos2] - v * allpass2->feedback2;
	v += allpass2->feedback1 * allpass2->buf1[allpass2->pos1];
	allpass2->buf2[allpass2->pos2] = allpass2->decay1 * allpass2->buf1[allpass2->pos1] -
		v * allpass2->feedback1;
	allpass2->buf1[allpass2->pos1] = v;
	allpass2->pos1 = (allpass2->pos1 + 1) % allpass2->size1;
	allpass2->pos2 = (allpass2->pos2 + 1) % allpass2->size2;
	return out;
}

static inline jfloat allpass2_get1(sf_rv_allpass2_st *allpass2, jint offset){
	if (offset > allpass2->size1)
		return allpass2->buf1[allpass2->pos1];
	else if (offset <= 0)
		offset = 1;
	jint rp = allpass2->pos1 - offset;
	if (rp < 0)
		rp += allpass2->size1;
	return allpass2->buf1[rp];
}

static inline jfloat allpass2_get2(sf_rv_allpass2_st *allpass2, jint offset){
	if (offset > allpass2->size2)
		return allpass2->buf2[allpass2->pos2];
	else if (offset <= 0)
		offset = 1;
	jint rp = allpass2->pos2 - offset;
	if (rp < 0)
		rp += allpass2->size2;
	return allpass2->buf2[rp];
}

//
// allpass3
//
static inline void allpass3_make(sf_rv_allpass3_st *allpass3, jint size1, jint msize1, jint size2,
	jint size3, jfloat feedback1, jfloat feedback2, jfloat feedback3, jfloat decay1, jfloat decay2,
	jfloat decay3){
	size1 = clampi(size1, 1, SF_REVERB_AP3S1);
	msize1 = clampi(msize1, 1, SF_REVERB_AP3M1);
	if (msize1 > size1)
		msize1 = size1;
	jint newsize = size1 + msize1;
	allpass3->rpos1 = (msize1 * 2) % newsize;
	allpass3->wpos1 = 0;
	allpass3->pos2 = 0;
	allpass3->pos3 = 0;
	allpass3->size1 = newsize;
	allpass3->msize1 = msize1;
	allpass3->size2 = clampi(size2, 1, SF_REVERB_AP3S2);
	allpass3->size3 = clampi(size3, 1, SF_REVERB_AP3S3);
	allpass3->feedback1 = feedback1;
	allpass3->feedback2 = feedback2;
	allpass3->feedback3 = feedback3;
	allpass3->decay1 = decay1;
	allpass3->decay2 = decay2;
	allpass3->decay3 = decay3;
	// memset(allpass3->buf1, 0, sizeof(jfloat) * allpass3->size1);
	// memset(allpass3->buf2, 0, sizeof(jfloat) * allpass3->size2);
	// memset(allpass3->buf3, 0, sizeof(jfloat) * allpass3->size3);
}

static inline jfloat allpass3_step(sf_rv_allpass3_st *allpass3, jfloat v, jfloat mod){
	mod = (mod + 1.0f) * (jfloat)allpass3->msize1;
	jfloat floormod = floorf(mod);
	jfloat mfrac = mod - floormod;
	jint rpos1 = allpass3->rpos1 - (jint)floormod;
	if (rpos1 < 0)
		rpos1 += allpass3->size1;
	jint rpos2 = rpos1 - 1;
	if (rpos2 < 0)
		rpos2 += allpass3->size1;
	v += allpass3->feedback3 * allpass3->buf3[allpass3->pos3];
	jfloat out = allpass3->decay3 * allpass3->buf3[allpass3->pos3] - allpass3->feedback3 * v;
	v += allpass3->feedback2 * allpass3->buf2[allpass3->pos2];
	allpass3->buf3[allpass3->pos3] = allpass3->decay2 * allpass3->buf2[allpass3->pos2] -
		allpass3->feedback2 * v;
	jfloat tmp = allpass3->buf1[rpos2] * mfrac + allpass3->buf1[rpos1] * (1.0f - mfrac);
	v += allpass3->feedback1 * tmp;
	allpass3->buf2[allpass3->pos2] = allpass3->decay1 * tmp - allpass3->feedback1 * v;
	allpass3->buf1[allpass3->wpos1] = v;
	allpass3->wpos1 = (allpass3->wpos1 + 1) % allpass3->size1;
	allpass3->rpos1 = (allpass3->rpos1 + 1) % allpass3->size1;
	allpass3->pos2 = (allpass3->pos2 + 1) % allpass3->size2;
	allpass3->pos3 = (allpass3->pos3 + 1) % allpass3->size3;
	return out;
}

static inline jfloat allpass3_get1(sf_rv_allpass3_st *allpass3, jint offset){
	if (offset > allpass3->size1)
		return allpass3->buf1[allpass3->rpos1];
	else if (offset <= 0)
		offset = 1;
	jint rp = allpass3->rpos1 - offset;
	if (rp < 0)
		rp += allpass3->size1;
	return allpass3->buf1[rp];
}

static inline jfloat allpass3_get2(sf_rv_allpass3_st *allpass3, jint offset){
	if (offset > allpass3->size2)
		return allpass3->buf2[allpass3->pos2];
	else if (offset <= 0)
		offset = 1;
	jint rp = allpass3->pos2 - offset;
	if (rp < 0)
		rp += allpass3->size2;
	return allpass3->buf2[rp];
}

static inline jfloat allpass3_get3(sf_rv_allpass3_st *allpass3, jint offset){
	if (offset > allpass3->size3)
		return allpass3->buf3[allpass3->pos3];
	else if (offset <= 0)
		offset = 1;
	jint rp = allpass3->pos3 - offset;
	if (rp < 0)
		rp += allpass3->size3;
	return allpass3->buf3[rp];
}

//
// allpassm
//
static inline void allpassm_make(sf_rv_allpassm_st *allpassm, jint size, jint msize, jfloat feedback,
	jfloat decay){
	size = clampi(size, 1, SF_REVERB_APMS);
	msize = clampi(msize, 1, SF_REVERB_APMM);
	if (msize > size)
		msize = size;
	jint newsize = size + msize;
	allpassm->rpos = (msize * 2) % newsize;
	allpassm->wpos = 0;
	allpassm->size = newsize;
	allpassm->msize = msize;
	allpassm->feedback = feedback;
	allpassm->decay = decay;
	allpassm->z1 = 0;
	// memset(allpassm->buf, 0, sizeof(jfloat) * allpassm->size);
}

static inline jfloat allpassm_step(sf_rv_allpassm_st *allpassm, jfloat v, jfloat mod, jfloat fbmod){
	jfloat mfeedback = allpassm->feedback + fbmod;
	mod = (mod + 1.0f) * (jfloat)allpassm->msize;
	jfloat floormod = floorf(mod);
	jfloat mfrac = 1.0f - mod + floormod;
	jint rpos1 = allpassm->rpos - (jint)floormod;
	if (rpos1 < 0)
		rpos1 += allpassm->size;
	jint rpos2 = rpos1 - 1;
	if (rpos2 < 0)
		rpos2 += allpassm->size;
	allpassm->z1 = allpassm->buf[rpos2] + mfrac * (allpassm->buf[rpos1] - allpassm->z1);
	allpassm->rpos = (allpassm->rpos + 1) % allpassm->size;
	allpassm->buf[allpassm->wpos] = v + allpassm->z1 * mfeedback;
	v = allpassm->decay * allpassm->z1 - allpassm->buf[allpassm->wpos] * mfeedback;
	allpassm->wpos = (allpassm->wpos + 1) % allpassm->size;
	return v;
}

//
// comb
//
static inline void comb_make(sf_rv_comb_st *comb, jint size){
	comb->pos = 0;
	comb->size = clampi(size, 1, SF_REVERB_CS);
	// memset(comb->buf, 0, sizeof(jfloat) * comb->size);
}

static inline jfloat comb_step(sf_rv_comb_st *comb, jfloat v, jfloat feedback){
	v = comb->buf[comb->pos] * feedback + v;
	comb->buf[comb->pos] = v;
	comb->pos = (comb->pos + 1) % comb->size;
	return v;
}

//
//
// reverb implementation
//
//

// now that all the components are done (thank god), we can start on the actual reverb effect

void sf_presetreverb(sf_reverb_state_st *rv, jint rate, sf_reverb_preset preset){
	// sorry for the bad formatting, I've tried to cram this in as best as I could
	struct {
		jint osf; jfloat p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15, p16;
	} ps[] = {

//OSF ERtoLt ERWet Dry ERFac ERWdth Wdth Wet Wander BassB Spin InpLP BasLP DmpLP OutLP RT60  Delay
{1, 0.40f, -9.0f,-10, 1.6f, 0.7f, 1.0f, -0, 0.27f, 0.15f, 0.7f,17000, 500, 7000,10000, 3.2f,0.020f},
{2, 0.30f, -9.0f, -8, 1.0f, 0.7f, 1.0f, -8, 0.32f, 0.25f, 0.7f,18000, 600, 9000,17000, 2.1f,0.010f},
{1, 0.30f, -9.0f, -8, 1.0f, 0.7f, 1.0f, -8, 0.27f, 0.20f, 0.5f,18000, 600, 7000, 9000, 2.3f,0.010f},
{2, 0.30f, -9.0f, -8, 1.2f, 0.7f, 1.0f, -8, 0.27f, 0.20f, 0.7f,18000, 500, 8000,16000, 2.8f,0.010f},
{1, 0.30f, -9.0f, -8, 1.2f, 0.7f, 1.0f, -8, 0.25f, 0.15f, 0.5f,18000, 500, 6000, 8000, 2.9f,0.010f},
{2, 0.20f, -9.0f, -8, 1.4f, 0.7f, 1.0f, -8, 0.17f, 0.20f, 1.0f,18000, 400, 9000,14000, 3.8f,0.018f},
{2, 0.20f, -9.0f, -8, 1.5f, 0.7f, 1.0f, -8, 0.20f, 0.20f, 0.5f,18000, 400, 5000, 7000, 4.2f,0.018f},
{2, 0.70f, -8.0f, -8, 0.7f,-0.4f, 0.8f, -8, 0.20f, 0.30f, 1.6f,18000,1000,18000,18000, 0.5f,0.005f},
{3, 0.70f, -8.0f, -8, 0.8f, 0.6f, 0.9f, -8, 0.30f, 0.30f, 0.4f,18000, 300,10000,18000, 0.5f,0.005f},
{2, 0.50f, -8.0f, -8, 1.2f,-0.4f, 0.8f, -8, 0.20f, 0.10f, 1.6f,18000,1000,18000,18000, 0.8f,0.008f},
{2, 0.50f, -8.0f, -8, 1.2f, 0.6f, 0.9f, -8, 0.30f, 0.10f, 0.4f,18000, 300,10000,18000, 1.2f,0.016f},
{2, 0.20f, -8.0f, -8, 2.2f,-0.4f, 0.9f, -8, 0.20f, 0.10f, 1.6f,18000,1000,16000,18000, 1.8f,0.010f},
{2, 0.20f, -8.0f, -8, 2.2f, 0.6f, 0.9f, -8, 0.30f, 0.10f, 0.4f,18000, 500, 9000,18000, 1.9f,0.020f},
{2, 0.50f, -7.0f, -7, 1.2f,-0.4f, 0.8f,-70, 0.20f, 0.10f, 1.6f,18000,1000,18000,18000, 0.8f,0.008f},
{2, 0.50f, -7.0f, -7, 1.2f, 0.6f, 0.9f,-70, 0.30f, 0.10f, 0.4f,18000, 300,10000,18000, 1.2f,0.016f},
{2, 0.00f,-70.0f,-20, 1.0f, 1.0f, 1.0f, -8, 0.20f, 0.10f, 1.6f,18000,1000,16000,18000, 1.8f,0.000f},
{2, 0.00f,-70.0f,-20, 1.0f, 1.0f, 1.0f, -8, 0.30f, 0.20f, 0.4f,18000, 500, 9000,18000, 1.9f,0.000f},
{2, 0.10f,-16.0f,-15, 1.0f, 0.1f, 1.0f, -5, 0.35f, 0.05f, 1.0f,18000, 100,10000,18000,12.0f,0.000f},
{2, 0.10f,-16.0f,-15, 1.0f, 0.1f, 1.0f, -5, 0.40f, 0.05f, 1.0f,18000, 100, 9000,18000,30.0f,0.000f}

	};

	#define CASE(prs, i)                                                                        \
		case prs: sf_advancereverb(rv, rate, ps[i].osf, ps[i].p1, ps[i].p2, ps[i].p3, ps[i].p4, \
			ps[i].p5, ps[i].p6, ps[i].p7, ps[i].p8, ps[i].p9, ps[i].p10, ps[i].p11, ps[i].p12,  \
			ps[i].p13, ps[i].p14, ps[i].p15, ps[i].p16); return;
	switch (preset){
		CASE(SF_REVERB_PRESET_DEFAULT    ,  0)
		CASE(SF_REVERB_PRESET_SMALLHALL1 ,  1)
		CASE(SF_REVERB_PRESET_SMALLHALL2 ,  2)
		CASE(SF_REVERB_PRESET_MEDIUMHALL1,  3)
		CASE(SF_REVERB_PRESET_MEDIUMHALL2,  4)
		CASE(SF_REVERB_PRESET_LARGEHALL1 ,  5)
		CASE(SF_REVERB_PRESET_LARGEHALL2 ,  6)
		CASE(SF_REVERB_PRESET_SMALLROOM1 ,  7)
		CASE(SF_REVERB_PRESET_SMALLROOM2 ,  8)
		CASE(SF_REVERB_PRESET_MEDIUMROOM1,  9)
		CASE(SF_REVERB_PRESET_MEDIUMROOM2, 10)
		CASE(SF_REVERB_PRESET_LARGEROOM1 , 11)
		CASE(SF_REVERB_PRESET_LARGEROOM2 , 12)
		CASE(SF_REVERB_PRESET_MEDIUMER1  , 13)
		CASE(SF_REVERB_PRESET_MEDIUMER2  , 14)
		CASE(SF_REVERB_PRESET_PLATEHIGH  , 15)
		CASE(SF_REVERB_PRESET_PLATELOW   , 16)
		CASE(SF_REVERB_PRESET_LONGREVERB1, 17)
		CASE(SF_REVERB_PRESET_LONGREVERB2, 18)
	}
	#undef CASE
}

void sf_advancereverb(sf_reverb_state_st *rv, jint rate,
	jint oversamplefactor, jfloat ertolate, jfloat erefwet, jfloat dry, jfloat ereffactor,
	jfloat erefwidth, jfloat width, jfloat wet, jfloat wander, jfloat bassb, jfloat spin, jfloat inputlpf,
	jfloat basslpf, jfloat damplpf, jfloat outputlpf, jfloat rt60, jfloat delay){

	rv->ertolate = ertolate;
	rv->erefwet = db2lin(erefwet);
	rv->dry = db2lin(dry);
	wet = db2lin(wet);
	rv->wet1 = wet * (width * 0.5f + 0.5f);
	rv->wet2 = wet * ((1.0f - width) * 0.5f);
	rv->wander = wander;
	rv->bassb = bassb;

	earlyref_make(&rv->earlyref, rate, ereffactor, erefwidth);

	oversample_make(&rv->oversampleL, oversamplefactor);
	rv->oversampleR = rv->oversampleL;
	jint osrate = rate * rv->oversampleL.factor;

	dccut_make(&rv->dccutL, osrate, 5.0f);
	rv->dccutR = rv->dccutL;

	noise_make(&rv->noise);

	lfo_make(&rv->lfo1, osrate, spin);
	iir1_makeLPF(&rv->lfo1_lpf, osrate, 20.0f);
	lfo_make(&rv->lfo2, osrate, sqrtf(100.0f - (10.0f - spin) * (10.0f - spin)) * 0.5f);
	iir1_makeLPF(&rv->lfo2_lpf, osrate, 12.0f);

	static const jint diffLc[10] = { 617, 535, 434, 347, 218, 162, 144, 122, 109, 74 };
	static const jint diffRc[10] = { 603, 547, 416, 364, 236, 162, 140, 131, 111, 79 };
	jint totfactor = osrate / 34125;
	jint msize = nextprime(10 * osrate / 34125);
	for (jint i = 0; i < 10; i++){
		allpassm_make(&rv->diffL[i], nextprime(diffLc[i] * totfactor), msize, -0.78f, 1);
		allpassm_make(&rv->diffR[i], nextprime(diffRc[i] * totfactor), msize, -0.78f, 1);
	}

	static const jint crossLc[4] = { 430, 341, 264, 174 };
	static const jint crossRc[4] = { 447, 324, 247, 191 };
	for (jint i = 0; i < 4; i++){
		allpass_make(&rv->crossL[i], nextprime(crossLc[i] * totfactor), 0.78f, 1);
		allpass_make(&rv->crossR[i], nextprime(crossRc[i] * totfactor), 0.78f, 1);
	}

	iir1_makeLPF(&rv->clpfL, osrate, inputlpf);
	rv->clpfR = rv->clpfL;

	delay_make(&rv->cdelayL , nextprime(1572 * totfactor));
	delay_make(&rv->cdelayR , nextprime(  16 * totfactor));
	delay_make(&rv->dampdL  , nextprime(   2 * totfactor));
	delay_make(&rv->dampdR  , nextprime(       totfactor));
	delay_make(&rv->cbassd1L, nextprime(1055 * totfactor));
	delay_make(&rv->cbassd1R, nextprime(1460 * totfactor));
	delay_make(&rv->cbassd2L, nextprime( 344 * totfactor));
	delay_make(&rv->cbassd2R, nextprime( 500 * totfactor));

	biquad_makeAPF(&rv->bassapL, osrate, 150.0f, 4.0f);
	rv->bassapR = rv->bassapL;

	biquad_makeLPF(&rv->basslpL, osrate, basslpf, 2.0f);
	rv->basslpR = rv->basslpL;

	iir1_makeLPF(&rv->damplpL, osrate, damplpf);
	rv->damplpR = rv->damplpL;

	jfloat decay0 = powf(10.0f, log10f(0.237f) / rt60);
	jfloat decay1 = powf(10.0f, log10f(0.938f) / rt60);
	jfloat decay2 = powf(10.0f, log10f(0.844f) / rt60);
	jfloat decay3 = powf(10.0f, log10f(0.906f) / rt60);
	rv->loopdecay = decay0;
	msize = nextprime(32 * totfactor);
	allpassm_make(&rv->dampap1L, nextprime(239 * totfactor), msize, 0.375f, decay2);
	allpassm_make(&rv->dampap1R, nextprime(205 * totfactor), msize, 0.375f, decay2);
	allpassm_make(&rv->dampap2L, nextprime(392 * totfactor), msize, 0.312f, decay3);
	allpassm_make(&rv->dampap2R, nextprime(329 * totfactor), msize, 0.312f, decay3);

	allpass2_make(&rv->cbassap1L, nextprime(1944 * totfactor), nextprime(612 * totfactor),
		0.250f, 0.406f, decay1, decay2);
	allpass2_make(&rv->cbassap1R, nextprime(2032 * totfactor), nextprime(368 * totfactor),
		0.250f, 0.406f, decay1, decay2);

	allpass3_make(&rv->cbassap2L,
		nextprime(1212 * totfactor),
		nextprime( 121 * totfactor),
		nextprime( 816 * totfactor),
		nextprime(1264 * totfactor),
		0.250f, 0.250f, 0.406f, decay1, decay1, decay2);
	allpass3_make(&rv->cbassap2R,
		nextprime(1452 * totfactor),
		nextprime(   5 * totfactor),
		nextprime( 688 * totfactor),
		nextprime(1340 * totfactor),
		0.250f, 0.250f, 0.406f, decay1, decay1, decay2);

	static const jint outco[32] = {
		  1,  40, 192, 276, 321, 110, 468, 1572, 121, 480, 103, 26, 780, 1200, 310, 780,
		625, 468, 312,  24,  36, 790, 189,    8,  10, 359,  30, 10, 109, 1310, 800,  10
	};
	for (jint i = 0; i < 32; i++)
		rv->outco[i] = outco[i] * totfactor;

	comb_make(&rv->combL, nextprime(22 * osrate / 1000));
	rv->combR = rv->combL;

	biquad_makeLPF(&rv->lastlpfL, osrate, outputlpf, 1.0f);
	rv->lastlpfR = rv->lastlpfL;

	jint delaysamp = osrate * delay;
	if (delaysamp >= 0){
		delay_make(&rv->inpdelayL, 0);
		delay_make(&rv->inpdelayR, 0);
		delay_make(&rv->lastdelayL, delaysamp);
		delay_make(&rv->lastdelayR, delaysamp);
	}
	else{
		delay_make(&rv->inpdelayL, -delaysamp);
		delay_make(&rv->inpdelayR, -delaysamp);
		delay_make(&rv->lastdelayL, 0);
		delay_make(&rv->lastdelayR, 0);
	}
}

void sf_reverb_process(sf_reverb_state_st *rv, jint size, sf_sample_st *input, sf_sample_st *output){
	// extra hardcoded constants
	const jfloat modnoise1 = 0.09f;
	const jfloat modnoise2 = 0.06f;
	const jfloat crossfeed = 0.4f;

	// oversample buffer
	jfloat osL[SF_REVERB_OF], osR[SF_REVERB_OF];

	for (jint i = 0; i < size; i++){
		// early reflection
		sf_sample_st er = earlyref_step(&rv->earlyref, input[i]);
		jfloat erL = er.L * rv->ertolate + input[i].L;
		jfloat erR = er.R * rv->ertolate + input[i].R;

		// oversample the single input jinto multiple outputs
		oversample_stepup(&rv->oversampleL, erL, osL);
		oversample_stepup(&rv->oversampleR, erR, osR);

		// for each oversampled sample...
		for (jint i2 = 0; i2 < rv->oversampleL.factor; i2++){
			// dc cut
			jfloat outL = dccut_step(&rv->dccutL, osL[i2]);
			jfloat outR = dccut_step(&rv->dccutR, osR[i2]);

			// noise
			jfloat mnoise = noise_step(&rv->noise);
			jfloat lfo = (lfo_step(&rv->lfo1) + modnoise1 * mnoise) * rv->wander;
			lfo = iir1_step(&rv->lfo1_lpf, lfo);
			mnoise *= modnoise2;

			// diffusion
			// for (jint i = 0, s = -1; i < 10; i++, s = -s){
			// 	outL = allpassm_step(&rv->diffL[i], outL, lfo * s, mnoise);
			// 	outR = allpassm_step(&rv->diffR[i], outR, lfo, mnoise * s);
			// }

			// cross fade
			jfloat crossL = outL, crossR = outR;
			for (jint i = 0; i < 4; i++){
				crossL = allpass_step(&rv->crossL[i], crossL);
				crossR = allpass_step(&rv->crossR[i], crossR);
			}
			outL = iir1_step(&rv->clpfL, outL + crossfeed * crossR);
			outR = iir1_step(&rv->clpfR, outR + crossfeed * crossL);

			// bass boost
			crossL = delay_getlast(&rv->cdelayL);
			crossR = delay_getlast(&rv->cdelayR);
			outL += rv->loopdecay *
				(crossR + rv->bassb * biquad_step(&rv->basslpL, biquad_step(&rv->bassapL, crossR)));
			outR += rv->loopdecay *
				(crossL + rv->bassb * biquad_step(&rv->basslpR, biquad_step(&rv->bassapR, crossL)));

			// dampening
			outL = allpassm_step(&rv->dampap2L,
				delay_step(&rv->dampdL,
				allpassm_step(&rv->dampap1L,
				iir1_step(&rv->damplpL, outL), lfo, mnoise)),
				-lfo, -mnoise);
			outR = allpassm_step(&rv->dampap2R,
				delay_step(&rv->dampdR,
				allpassm_step(&rv->dampap1R,
				iir1_step(&rv->damplpR, outR), -lfo, -mnoise)),
				lfo, mnoise);

			// update cross fade bass boost delay
			delay_step(&rv->cdelayL,
				allpass3_step(&rv->cbassap2L,
				delay_step(&rv->cbassd2L,
				allpass2_step(&rv->cbassap1L,
				delay_step(&rv->cbassd1L, outL))),
					lfo));
			delay_step(&rv->cdelayR,
				allpass3_step(&rv->cbassap2R,
				delay_step(&rv->cbassd2R,
				allpass2_step(&rv->cbassap1R,
				delay_step(&rv->cbassd1R, outR))),
					-lfo));

			//
			jfloat D1 =
				delay_get    (&rv->cbassd1L , rv->outco[ 0]);
			jfloat D2 =
				delay_get    (&rv->cbassd2L , rv->outco[ 1]) -
				delay_get    (&rv->cbassd2R , rv->outco[ 2]) +
				delay_get    (&rv->cbassd2L , rv->outco[ 3]) -
				delay_get    (&rv->cdelayR  , rv->outco[ 4]) -
				delay_get    (&rv->cbassd1R , rv->outco[ 5]) -
				delay_get    (&rv->cbassd2R , rv->outco[ 6]);
			jfloat D3 =
				delay_get    (&rv->cdelayL  , rv->outco[ 7]) +
				allpass2_get1(&rv->cbassap1L, rv->outco[ 8]) +
				allpass2_get2(&rv->cbassap1L, rv->outco[ 9]) -
				allpass2_get2(&rv->cbassap1R, rv->outco[10]) +
				allpass3_get1(&rv->cbassap2L, rv->outco[11]) +
				allpass3_get2(&rv->cbassap2L, rv->outco[12]) +
				allpass3_get3(&rv->cbassap2L, rv->outco[13]) -
				allpass3_get2(&rv->cbassap2R, rv->outco[14]);
			jfloat D4 =
				delay_get    (&rv->cdelayL  , rv->outco[15]);

			jfloat B1 =
				delay_get    (&rv->cbassd1R , rv->outco[16]);
			jfloat B2 =
				delay_get    (&rv->cbassd2R , rv->outco[17]) -
				delay_get    (&rv->cbassd2L , rv->outco[18]) +
				delay_get    (&rv->cbassd2R , rv->outco[19]) -
				delay_get    (&rv->cdelayL  , rv->outco[20]) -
				delay_get    (&rv->cbassd1L , rv->outco[21]) -
				delay_get    (&rv->cbassd2L , rv->outco[22]);
			jfloat B3 =
				delay_get    (&rv->cdelayR  , rv->outco[23]) +
				allpass2_get1(&rv->cbassap1R, rv->outco[24]) +
				allpass2_get2(&rv->cbassap1R, rv->outco[25]) -
				allpass2_get2(&rv->cbassap1L, rv->outco[26]) +
				allpass3_get1(&rv->cbassap2R, rv->outco[27]) +
				allpass3_get2(&rv->cbassap2R, rv->outco[28]) +
				allpass3_get3(&rv->cbassap2R, rv->outco[29]) -
				allpass3_get2(&rv->cbassap2L, rv->outco[30]);
			jfloat B4 =
				delay_get    (&rv->cdelayR  , rv->outco[31]);

			jfloat D = D1 * 0.469f + D2 * 0.219f + D3 * 0.064f + D4 * 0.045f;
			jfloat B = B1 * 0.469f + B2 * 0.219f + B3 * 0.064f + B4 * 0.045f;

			lfo = iir1_step(&rv->lfo2_lpf, lfo_step(&rv->lfo2) * rv->wander);
			outL = comb_step(&rv->combL, D, lfo);
			outR = comb_step(&rv->combR, B, -lfo);

			outL = delay_step(&rv->lastdelayL, biquad_step(&rv->lastlpfL, outL));
			outR = delay_step(&rv->lastdelayR, biquad_step(&rv->lastlpfR, outR));

			osL[i2] = outL * rv->wet1 + outR * rv->wet2 +
				delay_step(&rv->inpdelayL, osL[i2]) * rv->dry;
			osR[i2] = outR * rv->wet1 + outL * rv->wet2 +
				delay_step(&rv->inpdelayR, osR[i2]) * rv->dry;
		}

		jfloat outL = oversample_stepdown(&rv->oversampleL, osL);
		jfloat outR = oversample_stepdown(&rv->oversampleR, osR);
		outL += er.L * rv->erefwet + input[i].L * rv->dry;
		outR += er.R * rv->erefwet + input[i].R * rv->dry;
		output[i] = (sf_sample_st){ outL, outR };
	}
}