#### RemiAudio
#### Copyright (C) 2022-2024 Remilia Scarlet <remilia@posteo.jp>
#### Based on Zita-Rev1
#### Copyright (C) 2003-2017 Fons Adriaensen <fons@linuxaudio.org>
####
#### This program is free software; you can redistribute it and/or modify it
#### under the terms of the GNU General Public License as published by the Free
#### Software Foundation; either version 3 of the License, or (at your option)
#### any later version.
####
#### This program is distributed in the hope that it will be useful, but WITHOUT
#### ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
#### FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
#### more details.
####
#### You should have received a copy of the GNU General Public License along
#### with this program. If not, see <http://www.gnu.org/licenses/>.
require "./reverb"
require "./zita-paraeq"
require "./zita-presets.cr"
module RemiAudio::DSP
# Implements a Hall-like reverb effect.
#
# This is a port of Zita-Rev1, a mostly-a-hall-partially-a-plate reverb based
# on a feedback delay network and allpass comb filters.
#
# The original Zita-Rev1 supports an "ambisonic" mode, which allows for two
# inputs and four outputs. This implementation does not support the
# calculations for ambisonic mode, however, so it's strictly a reverb with
# stereo output.
class ZitaReverb < Reverb
PARA_EQ_MAX_CH = 2 # Change to 4 if we ever get ambisonic mode support
NUM_DIFFS = 8
NUM_FILTERS = 8
NUM_DELAYS = 8
ZITA_G = Math.sqrt(0.125)
DEFAULT_PRE_DELAY = 0.02
MINIMUM_PRE_DELAY = 0.02
MAXIMUM_PRE_DELAY = 0.10
DEFAULT_REVERB_XOVER = 369.0
MINIMUM_REVERB_XOVER = 50.0
MAXIMUM_REVERB_XOVER = 1000.0
DEFAULT_REVERB_TIME_LOW = 2.12
MINIMUM_REVERB_TIME_LOW = 1.0
MAXIMUM_REVERB_TIME_LOW = 8.0
DEFAULT_REVERB_TIME_HIGH = 2.84
MINIMUM_REVERB_TIME_HIGH = 1.0
MAXIMUM_REVERB_TIME_HIGH = 8.0
DEFAULT_REVERB_DAMPING = 4300.0
MINIMUM_REVERB_DAMPING = 1500.0
MAXIMUM_REVERB_DAMPING = 24000.0
DEFAULT_REVERB_EQ1_FREQ = 200.0
MINIMUM_REVERB_EQ1_FREQ = 40.0
MAXIMUM_REVERB_EQ1_FREQ = 2500.0
DEFAULT_REVERB_EQ1_GAIN = -1.3
MINIMUM_REVERB_EQ1_GAIN = -15.0
MAXIMUM_REVERB_EQ1_GAIN = 15.0
DEFAULT_REVERB_EQ2_FREQ = 2500.0
MINIMUM_REVERB_EQ2_FREQ = 160.0
MAXIMUM_REVERB_EQ2_FREQ = 10000.0
DEFAULT_REVERB_EQ2_GAIN = 0.0
MINIMUM_REVERB_EQ2_GAIN = -15.0
MAXIMUM_REVERB_EQ2_GAIN = 15.0
############################################################################
private class Diff
@i : Int32 = 0
@c : Float64
@line : Array(Float64)
def initialize(size : Int, @c : Float64)
@line = Array.new(size, 0.0f64)
end
@[AlwaysInline]
def process(x : Float64) : Float64
z = @line.unsafe_fetch(@i)
x -= @c * z
@line.unsafe_put(@i, x)
@i += 1
@i = 0 if @i == @line.size
z + @c * x
end
@[AlwaysInline]
def mute : Nil
@line.fill(0.0f64)
end
end
private class Filter
@gmf : Float64 = 0.0
@glo : Float64 = 0.0
@wlo : Float64 = 0.0
@whi : Float64 = 0.0
@slo : Float64 = 0.0
@shi : Float64 = 0.0
def initialize
end
@[AlwaysInline]
def process(x : Float64) : Float64
@slo += @wlo * (x - @slo) + 1.0e-10f64
x += @glo * @slo
@shi += @whi * (x - @shi)
@gmf * @shi
end
def setParams(del : Float64, tmf : Float64, tlo : Float64, newWlo : Float64, thi : Float64, chi : Float64) : Nil
@gmf = 0.001f64 ** (del / tmf)
@glo = (0.001f64 ** (del / tlo)) / @gmf - 1.0
@wlo = newWlo
g : Float64 = (0.001f64 ** (del / thi)) / @gmf
i : Float64 = (1 - g * g) / (2 * g * g * chi)
@whi = (Math.sqrt(1 + 4 * i) - 1.0) / (i + i)
end
def reset
@gmf = 0.0
@glo = 0.0
@wlo = 0.0
@whi = 0.0
@slo = 0.0
@shi = 0.0
end
end
private class Delay
@i : Int32 = 0
@line : Array(Float64)
def initialize(size : Int)
@line = Array.new(size, 0.0f64)
end
@[AlwaysInline]
def read : Float64
@line.unsafe_fetch(@i)
end
@[AlwaysInline]
def write(x : Float64) : Nil
@line.unsafe_put(@i, x)
@i += 1
@i = 0 if @i == @line.size
end
@[AlwaysInline]
def mute : Nil
@line.fill(0.0f64)
end
end
private class VDelay
@ir : Int32 = 0
@iw : Int32 = 0
@line : Array(Float64)
def initialize(size : Int)
@line = Array.new(size, 0.0f64)
end
@[AlwaysInline]
def read : Float64
ret = @line.unsafe_fetch(@ir)
@ir += 1
@ir = 0 if @ir == @line.size
ret
end
@[AlwaysInline]
def write(x : Float64) : Nil
@line.unsafe_put(@iw, x)
@iw += 1
@iw = 0 if @iw == @line.size
end
@[AlwaysInline]
def setDelay(del : Int32) : Nil
@ir = @iw - del
@ir += @line.size if @ir < 0
end
@[AlwaysInline]
def mute : Nil
@line.fill(0.0f64)
end
end
############################################################################
@sampleRate : Float64 = 44100.0
@invSampleRate : Float64 = 0.0
#@ambisonicMode : Bool = false
@vdelay0 : VDelay
@vdelay1 : VDelay
@diff1 : Array(Diff) = [] of Diff
@filt1 : Array(Filter) = [Filter.new, Filter.new, Filter.new, Filter.new,
Filter.new, Filter.new, Filter.new, Filter.new]
@delay : Array(Delay) = [] of Delay
@cntA1 : Int32 = 1
@cntB1 : Int32 = 1
@cntC1 : Int32 = 1
@cntA2 : Int32 = 0
@cntB2 : Int32 = 0
@cntC2 : Int32 = 0
@ipdel : Float64 = DEFAULT_PRE_DELAY
@xover : Float64 = DEFAULT_REVERB_XOVER
@rtlow : Float64 = DEFAULT_REVERB_TIME_LOW
@rtmid : Float64 = DEFAULT_REVERB_TIME_HIGH
@damping : Float64 = DEFAULT_REVERB_DAMPING
@opmix : Float64 = 1.0 # The "mix" value. Keep this at 1.0 since we only use this as an insert.
#@rgxyz : Float64 = 0.0 # For ambisonics, not used
@g0 : Float64 = 0.0
@d0 : Float64 = 0.0
@g1 : Float64 = 0.0
@d1 : Float64 = 0.0
@eq1 : ParametricEQ = ParametricEQ.new
@eq2 : ParametricEQ = ParametricEQ.new
TDIFF = [20346.0e-6f64,
24421.0e-6f64,
31604.0e-6f64,
27333.0e-6f64,
22904.0e-6f64,
29291.0e-6f64,
13458.0e-6f64,
19123.0e-6f64]
TDELAY = [153129.0e-6f64,
210389.0e-6f64,
127837.0e-6f64,
256891.0e-6f64,
174713.0e-6f64,
192303.0e-6f64,
125000.0e-6f64,
219991.0e-6f64]
def initialize(newSampleRate : Int32, numFrames : Int32)
@sampleRate = newSampleRate.to_f64
@invSampleRate = 1.0f64 / @sampleRate
# Create pre-delays
@vdelay0 = VDelay.new((0.1f64 * @sampleRate).to_i32!)
@vdelay1 = VDelay.new((0.1f64 * @sampleRate).to_i32!)
# Create main delay network
@diff1 = [] of Diff
@delay = [] of Delay
NUM_DIFFS.times do |i|
k1 = (TDIFF[i] * @sampleRate + 0.5).floor.to_i64!
k2 = (TDELAY[i] * @sampleRate + 0.5).floor.to_i64!
@diff1 << Diff.new(k1, ((i & 1) != 0 ? -0.6f64 : 0.6f64))
@delay << Delay.new(k2 - k1)
end
# Setup parametric EQs
@eq1.sampleRate = @sampleRate
@eq2.sampleRate = @sampleRate
# Prepare internal values
prepare(numFrames)
@eq1.prepare(numFrames)
@eq2.prepare(numFrames)
end
private def prepare(numFrames : Int32) : Nil
a : Int32 = @cntA1
b : Int32 = @cntB1
c : Int32 = @cntC1
@d0 = 0.0
@d1 = 0.0
k : Int32 = (((@ipdel - 0.02) * @sampleRate) + 0.5).trunc.to_i32!
@vdelay0.setDelay(k)
@vdelay1.setDelay(k)
@cntA2 = a
wlo : Float64 = 6.2832 * @xover * @invSampleRate
chi : Float64 = if @damping > (0.49 * @sampleRate)
2.0
else
1 - RemiMath.fastCos(6.2832 * @damping * @invSampleRate)
end
halfRtMid : Float64 = 0.5 * @rtmid
8.times do |i|
@filt1.unsafe_fetch(i).setParams(TDELAY.unsafe_fetch(i), @rtmid, @rtlow, wlo, halfRtMid, chi)
end
@cntB2 = b
t1 : Float64 = 0.7 / Math.sqrt(@rtmid)
@d0 = (-@g0) / numFrames
@d1 = (t1 - @g1) / numFrames
@cntC2 = c
end
def damping : Float64
@damping
end
def damping=(value : Float64)
ZitaReverb.checkDamping(value)
@damping = value
@cntB1 += 1
end
def predelay : Float64
@ipdel
end
def predelay=(value : Float64)
ZitaReverb.checkPredelay(value)
@ipdel = value
@cntA1 += 1
end
def crossover : Float64
@xover
end
def crossover=(value : Float64)
ZitaReverb.checkCrossover(value)
@xover = value
@cntB1 += 1
end
def timeLow : Float64
@rtlow
end
def timeLow=(value : Float64)
ZitaReverb.checkTimeLow(value)
@rtlow = value
@cntB1 += 1
end
def timeHigh : Float64
@rtmid
end
def timeHigh=(value : Float64)
ZitaReverb.checkTimeHigh(value)
@rtmid = value
@cntB1 += 1
@cntC1 += 1
end
def eq(number : Int) : Tuple(Float64, Float64)
if number == 1
{@eq1.f0, @eq1.g0}
elsif number == 2
{@eq2.f0, @eq2.g0}
else
raise "Bad EQ number: #{number}"
end
end
def setEq(number : Int, frequency : Float64, gain : Float64)
case number
when 1
ZitaReverb.checkEQ1Freq(frequency)
ZitaReverb.checkEQ1Gain(gain)
@eq1.setParam(frequency, gain)
when 2
ZitaReverb.checkEQ2Freq(frequency)
ZitaReverb.checkEQ2Gain(gain)
@eq2.setParam(frequency, gain)
else
raise "Bad EQ number: #{number}"
end
end
# Checks if `value` is a valid predelay time. If it is, this does nothing.
# Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkPredelay(value) : Nil
unless value >= MINIMUM_PRE_DELAY && value <= MAXIMUM_PRE_DELAY
raise ReverbError.new("Pre-delay for Zita must be between #{MINIMUM_PRE_DELAY} and #{MAXIMUM_PRE_DELAY}")
end
end
# Checks if `value` is a valid crossover frequency for the two bands. If it
# is, this does nothing. Otherwise it raises a `ReverbError` with an
# explanation.
@[AlwaysInline]
def self.checkCrossover(value) : Nil
unless value >= MINIMUM_REVERB_XOVER && value <= MAXIMUM_REVERB_XOVER
raise ReverbError.new("Crossover frequency for Zita must be between #{MINIMUM_REVERB_XOVER} and #{MAXIMUM_REVERB_XOVER}")
end
end
# Checks if `value` is a valid time length for the low band. If it is, this
# does nothing. Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkTimeLow(value) : Nil
unless value >= MINIMUM_REVERB_TIME_LOW && value <= MAXIMUM_REVERB_TIME_LOW
raise ReverbError.new("Low band time for Zita must be between #{MINIMUM_REVERB_TIME_LOW} and #{MAXIMUM_REVERB_TIME_LOW}")
end
end
# Checks if `value` is a valid time length for the high band. If it is,
# this does nothing. Otherwise it raises a `ReverbError` with an
# explanation.
@[AlwaysInline]
def self.checkTimeHigh(value) : Nil
unless value >= MINIMUM_REVERB_TIME_HIGH && value <= MAXIMUM_REVERB_TIME_HIGH
raise ReverbError.new("High band time for Zita must be between #{MINIMUM_REVERB_TIME_HIGH} and #{MAXIMUM_REVERB_TIME_HIGH}")
end
end
# Checks if `value` is a valid damping amount. If it is, this does nothing.
# Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkDamping(value) : Nil
unless value >= MINIMUM_REVERB_DAMPING && value <= MAXIMUM_REVERB_DAMPING
raise ReverbError.new("Damping for Zita must be between #{MINIMUM_REVERB_DAMPING} and #{MAXIMUM_REVERB_DAMPING}")
end
end
# Checks if `value` is a valid frequency for EQ1. If it is, this does nothing.
# Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkEQ1Freq(value) : Nil
unless value >= MINIMUM_REVERB_EQ1_FREQ && value <= MAXIMUM_REVERB_EQ1_FREQ
raise ReverbError.new("EQ1 frequency for Zita must be between #{MINIMUM_REVERB_EQ1_FREQ} and #{MAXIMUM_REVERB_EQ1_FREQ}")
end
end
# Checks if `value` is a valid gain for EQ1. If it is, this does nothing.
# Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkEQ1Gain(value) : Nil
unless value >= MINIMUM_REVERB_EQ1_GAIN && value <= MAXIMUM_REVERB_EQ1_GAIN
raise ReverbError.new("EQ1 gain for Zita must be between #{MINIMUM_REVERB_EQ1_GAIN} and #{MAXIMUM_REVERB_EQ1_GAIN}")
end
end
# Checks if `value` is a valid frequency for EQ2. If it is, this does nothing.
# Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkEQ2Freq(value) : Nil
unless value >= MINIMUM_REVERB_EQ2_FREQ && value <= MAXIMUM_REVERB_EQ2_FREQ
raise ReverbError.new("EQ2 frequency for Zita must be between #{MINIMUM_REVERB_EQ2_FREQ} and #{MAXIMUM_REVERB_EQ2_FREQ}")
end
end
# Checks if `value` is a valid gain for EQ2. If it is, this does nothing.
# Otherwise it raises a `ReverbError` with an explanation.
@[AlwaysInline]
def self.checkEQ2Gain(value) : Nil
unless value >= MINIMUM_REVERB_EQ2_GAIN && value <= MAXIMUM_REVERB_EQ2_GAIN
raise ReverbError.new("EQ2 gain for Zita must be between #{MINIMUM_REVERB_EQ2_GAIN} and #{MAXIMUM_REVERB_EQ2_GAIN}")
end
end
def mute : Nil
@vdelay0.mute
@vdelay1.mute
@diff1.each &.mute
@delay.each &.mute
@filt1.each &.reset
end
macro processZ(z, x0, x1)
{{z}} = {{x0}} - {{x1}}
{{x0}} += {{x1}}
{{x1}} = {{z}}
end
# Purposely unhygenic. Sets up variables for use in the process() methods.
private macro processCommon
z : Float64 = 0.0
x0 : Float64 = 0.0
x1 : Float64 = 0.0
x2 : Float64 = 0.0
x3 : Float64 = 0.0
x4 : Float64 = 0.0
x5 : Float64 = 0.0
x6 : Float64 = 0.0
x7 : Float64 = 0.0
end
private macro processLeftRight(left, right, outLeft, outRight)
# Write to the pre-delays
@vdelay0.write({{left}})
@vdelay1.write({{right}})
# Pre-Delay -> Diffuser
z = 0.375f64 * @vdelay0.read # Remi: This was 0.3 in the original Zita-Rev1 code
x0 = @diff1.unsafe_fetch(0).process(@delay.unsafe_fetch(0).read + z)
x1 = @diff1.unsafe_fetch(1).process(@delay.unsafe_fetch(1).read + z)
x2 = @diff1.unsafe_fetch(2).process(@delay.unsafe_fetch(2).read - z)
x3 = @diff1.unsafe_fetch(3).process(@delay.unsafe_fetch(3).read - z)
z = 0.375f64 * @vdelay1.read # Remi: This was 0.3 in the original Zita-Rev1 code
x4 = @diff1.unsafe_fetch(4).process(@delay.unsafe_fetch(4).read + z)
x5 = @diff1.unsafe_fetch(5).process(@delay.unsafe_fetch(5).read + z)
x6 = @diff1.unsafe_fetch(6).process(@delay.unsafe_fetch(6).read - z)
x7 = @diff1.unsafe_fetch(7).process(@delay.unsafe_fetch(7).read - z)
# Send signals through delay lines
processZ(z, x0, x1)
processZ(z, x2, x3)
processZ(z, x4, x5)
processZ(z, x6, x7)
processZ(z, x0, x2)
processZ(z, x1, x3)
processZ(z, x4, x6)
processZ(z, x5, x7)
processZ(z, x0, x4)
processZ(z, x1, x5)
processZ(z, x2, x6)
processZ(z, x3, x7)
# Mix
@g1 += @d1
{{outLeft}} = @g1 * (x1 + x2)
{{outRight}} = @g1 * (x1 - x2)
# Update FDN
@delay.unsafe_fetch(0).write(@filt1.unsafe_fetch(0).process(ZITA_G * x0))
@delay.unsafe_fetch(1).write(@filt1.unsafe_fetch(1).process(ZITA_G * x1))
@delay.unsafe_fetch(2).write(@filt1.unsafe_fetch(2).process(ZITA_G * x2))
@delay.unsafe_fetch(3).write(@filt1.unsafe_fetch(3).process(ZITA_G * x3))
@delay.unsafe_fetch(4).write(@filt1.unsafe_fetch(4).process(ZITA_G * x4))
@delay.unsafe_fetch(5).write(@filt1.unsafe_fetch(5).process(ZITA_G * x5))
@delay.unsafe_fetch(6).write(@filt1.unsafe_fetch(6).process(ZITA_G * x6))
@delay.unsafe_fetch(7).write(@filt1.unsafe_fetch(7).process(ZITA_G * x7))
end
def process(inputLeft : Array(Float64)|Slice(Float64), inputRight : Array(Float64)|Slice(Float64),
outputLeft : Array(Float64)|Slice(Float64), outputRight : Array(Float64)|Slice(Float64)) : Nil
unless inputLeft.size == outputLeft.size &&
inputRight.size == outputRight.size &&
inputLeft.size == inputRight.size
raise BufferSizeMismatchError.new("The length of the buffers must match")
end
prepare(inputLeft.size)
processCommon
inputLeft.size.times do |i|
processLeftRight(inputLeft[i], inputRight[i], outputLeft[i], outputRight[i])
end
# Run parametric EQs on the output
@eq1.process(outputLeft, outputRight)
@eq2.process(outputLeft, outputRight)
end
def process(inputLeft : Array(Float32)|Slice(Float32), inputRight : Array(Float32)|Slice(Float32),
outputLeft : Array(Float32)|Slice(Float32), outputRight : Array(Float32)|Slice(Float32)) : Nil
unless inputLeft.size == outputLeft.size &&
inputRight.size == outputRight.size &&
inputLeft.size == inputRight.size
raise BufferSizeMismatchError.new("The length of the buffers must match")
end
prepare(inputLeft.size)
processCommon
outLeft : Float64 = 0.0
outRight : Float64 = 0.0
inputLeft.size.times do |i|
processLeftRight(inputLeft[i].to_f64!, inputRight[i].to_f64!, outLeft, outRight)
outputLeft[i] = outLeft.to_f32!
outputRight[i] = outRight.to_f32!
end
# Run parametric EQs on the output
@eq1.process(outputLeft, outputRight)
@eq2.process(outputLeft, outputRight)
end
def process(input : Array(Float64)|Slice(Float64), output : Array(Float64)|Slice(Float64)) : Nil
unless input.size == output.size
raise BufferSizeMismatchError.new("The length of the buffers must match")
end
prepare(input.size.tdiv(2))
processCommon
outLeft : Float64 = 0.0
outRight : Float64 = 0.0
i = 0
while i < input.size
processLeftRight(input[i], input[i + 1], outLeft, outRight)
output[i] = outLeft * amount
output[i + 1] = outRight * amount
i += 2
end
# Run parametric EQs on the output
@eq1.process(output)
@eq2.process(output)
end
def process(input : Array(Float32)|Slice(Float32), output : Array(Float32)|Slice(Float32)) : Nil
unless input.size == output.size
raise BufferSizeMismatchError.new("The length of the buffers must match")
end
prepare(input.size.tdiv(2))
processCommon
outLeft : Float64 = 0.0
outRight : Float64 = 0.0
i = 0
while i < input.size
processLeftRight(input[i].to_f64!, input[i + 1].to_f64!, outLeft, outRight)
output[i] = (outLeft * amount).to_f32!
output[i + 1] = (outRight * amount).to_f32!
i += 2
end
# Run parametric EQs on the output
@eq1.process(output)
@eq2.process(output)
end
def process(buffer : Array(Float64)|Slice(Float64), amount : Float64|Float64 = 1.0) : Nil
prepare(buffer.size.tdiv(2))
processCommon
outLeft : Float64 = 0.0
outRight : Float64 = 0.0
i = 0
while i < buffer.size
processLeftRight(buffer[i], buffer[i + 1], outLeft, outRight)
buffer[i] += (outLeft * amount)
buffer[i + 1] += (outRight * amount)
i += 2
end
# Run parametric EQs on the output
@eq1.process(buffer)
@eq2.process(buffer)
end
def process(buffer : Array(Float32)|Slice(Float32), amount : Float32|Float64 = 1.0) : Nil
prepare(buffer.size.tdiv(2))
processCommon
outLeft : Float64 = 0.0
outRight : Float64 = 0.0
i = 0
while i < buffer.size
processLeftRight(buffer[i].to_f64!, buffer[i + 1].to_f64!, outLeft, outRight)
buffer[i] += (outLeft * amount).to_f32!
buffer[i + 1] += (outRight * amount).to_f32!
i += 2
end
# Run parametric EQs on the output
@eq1.process(buffer)
@eq2.process(buffer)
end
def updateInternalParams(*, numFrames : Int32 = 0)
prepare(numFrames)
end
def usePreset(preset : Reverb::Preset) : Nil
if data = preset.as?(ZitaReverb::Preset)
self.damping = data.damping
self.predelay = data.predelay
self.crossover = data.crossover
self.timeLow = data.timeLow
self.timeHigh = data.timeHigh
self.setEq(1, data.eq1Freq, data.eq1Gain)
self.setEq(2, data.eq2Freq, data.eq2Gain)
else
raise "Bad preset type for Zita: #{preset.class}"
end
end
end
end