/*
m_sg_reverb_q - A quadraphonic reverb based on SCReverb

DESCRIPTION
This is a four channel reverb based on the original SCReverb code by Sean
Costello.

SYNTAX
aOut1, aOut2, aOut3, aOut4 m_sg_reverb_q aIn1, aIn2, aIn3, aIn4, iPredelay, \
	kFeedback, kCutoff, kMod[, kLowCut, kMix]

INITIALIZATION
iPredelay - Predelay of the reverb in seconds. Internally limited to the range
	of 0 to .2

PERFORMANCE
aOut1, aOut2, aOut3, aOut4 - The four audio outputs
aIn1, aIn2, aIn3, aIn4 - The four audio inputs
kFeedback - The feedback level of the delay lines, i.e. how long it should
	reverberate, between 0 and 1
kCutoff - High frequency damping, the cutoff frequency of a simple 1st order
	lowpass filter, between 100 and 20000. The lower the cutoff the quicker
	the reverb will decay
kMod - Modulation of the delays, between 0 and 10. Usually, 1 is a good
	starting point
kLowCut (optional) - Optional cutoff frequency for low frequency rolloff.
	The default of -1 will roll off at 20Hz
	(Default: -1)
kMix (optional) - Mix level of the reverb, where 0 is completely dry and 1
	completely wet. This default to 1
	(Default: 1)

CREDITS
Author: Jeanette C., based on the code by Sean Costello
*/

opcode m_sg_reverb_q, aaaa, aaaaikkkJP
	aIn1, aIn2, aIn3, aIn4, iPredelay, kFeedback, kCutoff, kMod, \
	kLowCut, kMix xin

	; Setup and initialisation
	; Limit the input parameters
	iPredelay = limit:i(iPredelay, 0, .2)
	kFeedback = limit:k(kFeedback, 0, 1)
	kCutoff = limit:k(kCutoff, 100, 20000)
	kMod = limit:k(kMod, 0, 10)
	if (kLowCut <0) then
		kLowCut = 20
	else:
		kLowCut = limit:k(kLowCut, 0, (sr*.45))
	endif
	kMix = limit:k(kMix, 0, 1)
	iInVol init .7
	iOutVol init .7

	; Initialise the output signals
	aOut1 init 0
	aOut2 init 0
	aOut3 init 0
	aOut4 init 0

	; The four predelayed versions
	aPre1 init 0
	aPre2 init 0
	aPre3 init 0
	aPre4 init 0

	; Denormalise the inputs (just in case)
	denorm(aIn1, aIn2, aIn3, aIn4)

	; Initialise the delayline signals
	aFilt[] init 16 ; 16 delaylines
	kRand[] init 16 ; a randomiser per delayline

	; Set up the initial length of the delay lines
	if (sr == 48000) then
		iDelTimes[] = fillarray(2689, 3011, 3499, 3863, 4243, 4483, 2311, 2099, \
			3259, 2213, 2851, 2503, 3677, 4373, 1777, 4057)
	elseif (sr == 88200) then
		iDelTimes[] = fillarray(4943, 5531, 6427, 7109, 7793, 8243, 4283, 3863, \
			5987, 4073, 5237, 4603, 6763, 8017, 3271, 7459)
	elseif (sr == 96000) then
		iDelTimes[] = fillarray(5381, 6011, 7001, 7741, 8501, 8971, 4663, 4201, \
			6521, 4423, 5701, 5023, 7369, 8747, 3559, 8123)
	else ; 44.1 and everything else
		iDelTimes[] = fillarray(2473, 2767, 3217, 3557, 3907, 4127, 2143, 1933, \
			2999, 2039, 2621, 2309, 3389, 4019, 1637, 3733)
	endif

	; Rescale the delay times
	iDelTimes /= sr

	; Add the predelay
	if (iPredelay == 0) then
		aPre1 = aIn1
		aPre2 = aIn2
		aPre3 = aIn3
		aPre4 = aIn4
	else
		aPre1 = delay(aIn1, iPredelay)
		aPre2 = delay(aIn2, iPredelay)
		aPre3 = delay(aIn3, iPredelay)
		aPre4 = delay(aIn4, iPredelay)
	endif
	aPre1 *= iInVol
	aPre2 *= iInVol
	aPre3 *= iInVol
	aPre4 *= iInVol

	; Set up the randomisers
	kRand[0] = randi:k(.0013, 3.109, .06)
	kRand[1] = randi:k(.0011, 3.511, .9)
	kRand[2] = randi:k(.0017, 1.109, .7)
	kRand[3] = randi:k(.0007, 3.967, .3)
	kRand[4] = randi:k(.0013, 2.341, .63)
	kRand[5] = randi:k(.0011, 1.901, .7)
	kRand[6] = randi:k(.0017, 0.887, .9)
	kRand[7] = randi:k(.0007, 3.221, .44)
	kRand[8] = randi:k(.0005, .823, .03)
	kRand[9] = randi:k(.0029, 2.393, .47)
	kRand[10] = randi:k(.0019, 3.334, .54)
	kRand[11] = randi:k(.0023, 2.729, .28)
	kRand[12] = randi:k(.0005, 3.919, .8)
	kRand[13] = randi:k(.0029, 3.463, .73)
	kRand[14] = randi:k(.0017, 2.857, .36)
	kRand[15] = randi:k(.0023, 3.623, .67)

	; Calculate the junction pressure of the 16 delaylines
	apj = .125 * (aFilt[0] + aFilt[1] + aFilt[2] + aFilt[3] + aFilt[4] \
		+ aFilt[5] + aFilt[6] + aFilt[7] + aFilt[8] + aFilt[9] + aFilt[10] \
		+ aFilt[11] + aFilt[12] + aFilt[13] + aFilt[14] + aFilt[15])
	
	; Setup the real delay network
	aDummy[] init 16
	aDelays[] init 16
	aDummy[0] delayr 1
	aDelays[0] = deltapi(iDelTimes[0] + (kRand[0] * kMod))
	delayw(aIn1 + apj - aFilt[0])
	aDummy[1] delayr 1
	aDelays[1] = deltapi(iDelTimes[1] + (kRand[1] * kMod))
	delayw(aIn2 + apj - aFilt[1])
	aDummy[2] delayr 1
	aDelays[2] = deltapi(iDelTimes[2] + (kRand[2] * kMod))
	delayw(aIn1 + apj - aFilt[2] - aFilt[10])
	aDummy[3] delayr 1
	aDelays[3] = deltapi(iDelTimes[3] + (kRand[3] * kMod))
	delayw(aIn2 + apj - aFilt[3] - aFilt[11])
	aDummy[4] delayr 1
	aDelays[4] = deltapi(iDelTimes[4] + (kRand[4] * kMod))
	delayw(aIn1 + apj - aFilt[4])
	aDummy[5] delayr 1
	aDelays[5] = deltapi(iDelTimes[5] + (kRand[5] * kMod))
	delayw(aIn2 + apj - aFilt[5])
	aDummy[6] delayr 1
	aDelays[6] = deltapi(iDelTimes[6] + (kRand[6] * kMod))
	delayw(aIn1 + apj - aFilt[6] - aFilt[14])
	aDummy[7] delayr 1
	aDelays[7] = deltapi(iDelTimes[7] + (kRand[7] * kMod))
	delayw(aIn2 + apj - aFilt[7] - aFilt[15])
	aDummy[8] delayr 1
	aDelays[8] = deltapi(iDelTimes[8] + (kRand[8] * kMod))
	delayw(aIn3 + apj - aFilt[8] - aFilt[0])
	aDummy[9] delayr 1
	aDelays[9] = deltapi(iDelTimes[9] + (kRand[9] * kMod))
	delayw(aIn4 + apj - aFilt[9] - aFilt[1])
	aDummy[10] delayr 1
	aDelays[10] = deltapi(iDelTimes[10] + (kRand[10] * kMod))
	delayw(aIn3 + apj - aFilt[10])
	aDummy[11] delayr 1
	aDelays[11] = deltapi(iDelTimes[11] + (kRand[11] * kMod))
	delayw(aIn4 + apj - aFilt[11])
	aDummy[12] delayr 1
	aDelays[12] = deltapi(iDelTimes[12] + (kRand[12] * kMod))
	delayw(aIn3 + apj - aFilt[12] - aFilt[4])
	aDummy[13] delayr 1
	aDelays[13] = deltapi(iDelTimes[13] + (kRand[13] * kMod))
	delayw(aIn4 + apj - aFilt[13] - aFilt[5])
	aDummy[14] delayr 1
	aDelays[14] = deltapi(iDelTimes[14] + (kRand[14] * kMod))
	delayw(aIn3 + apj - aFilt[14])
	aDummy[15] delayr 1
	aDelays[15] = deltapi(iDelTimes[15] + (kRand[15] * kMod))
	delayw(aIn4 + apj - aFilt[15])

	; Apply the damping
	aFilt[0] = atone(tone((aDelays[0] * kFeedback), kCutoff), kLowCut)
	aFilt[1] = atone(tone((aDelays[1] * kFeedback), kCutoff), kLowCut)
	aFilt[2] = atone(tone((aDelays[2] * kFeedback), kCutoff), kLowCut)
	aFilt[3] = atone(tone((aDelays[3] * kFeedback), kCutoff), kLowCut)
	aFilt[4] = atone(tone((aDelays[4] * kFeedback), kCutoff), kLowCut)
	aFilt[5] = atone(tone((aDelays[5] * kFeedback), kCutoff), kLowCut)
	aFilt[6] = atone(tone((aDelays[6] * kFeedback), kCutoff), kLowCut)
	aFilt[7] = atone(tone((aDelays[7] * kFeedback), kCutoff), kLowCut)
	aFilt[8] = atone(tone((aDelays[8] * kFeedback), kCutoff), kLowCut)
	aFilt[9] = atone(tone((aDelays[9] * kFeedback), kCutoff), kLowCut)
	aFilt[10] = atone(tone((aDelays[10] * kFeedback), kCutoff), kLowCut)
	aFilt[11] = atone(tone((aDelays[11] * kFeedback), kCutoff), kLowCut)
	aFilt[12] = atone(tone((aDelays[12] * kFeedback), kCutoff), kLowCut)
	aFilt[13] = atone(tone((aDelays[13] * kFeedback), kCutoff), kLowCut)
	aFilt[14] = atone(tone((aDelays[14] * kFeedback), kCutoff), kLowCut)
	aFilt[15] = atone(tone((aDelays[15] * kFeedback), kCutoff), kLowCut)

	; Sum the delays for output
	aWet1 = dcblock2(aFilt[0] + aFilt[2] + aFilt[4] + aFilt[6]) * iOutVol
	aWet2 = dcblock2(aFilt[1] + aFilt[3] + aFilt[5] + aFilt[7]) * iOutVol
	aWet3 = dcblock2(aFilt[8] + aFilt[10] + aFilt[12] + aFilt[14]) * iOutVol
	aWet4 = dcblock2(aFilt[9] + aFilt[11] + aFilt[13] + aFilt[15]) * iOutVol
	aOut1 = ntrpol:a(aIn1, aWet1, kMix)
	aOut2 = ntrpol:a(aIn2, aWet2, kMix)
	aOut3 = ntrpol:a(aIn3, aWet3, kMix)
	aOut4 = ntrpol:a(aIn4, aWet4, kMix)

	; Output the audio
	xout(aOut1, aOut2, aOut3, aOut4)
endop

Csound mailing list
Csound@listserv.heanet.ie
https://listserv.heanet.ie/cgi-bin/wa?A0=CSOUND
Send bugs reports to
        https://github.com/csound/csound/issues
Discussions of bugs and features can be posted here