audio/ambivalence
2
0
Fork 0
mirror of https://codeberg.org/armin/ambivalence.git synced 2026-09-01 04:10:48 +02:00
Code Issues Projects Releases Packages Wiki Activity

init

Browse source
This commit is contained in:
Armin 2026-08-15 15:36:28 +02:00
commit 42cf8432bf
57 changed files with 5782 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

168
Source/DSP/AcousticMetrics.cpp Normal file
View file

@ -0,0 +1,168 @@
#include "AcousticMetrics.h"
#include <algorithm>
#include <cmath>
namespace FDNReverb {
void AcousticMetrics::prepare(double sr, float windowMs) {
sampleRate = sr;
analysisWindowMs = windowMs;
// sample rate time sample count
samples50ms = static_cast<int>(0.050 * sr);
samples80ms = static_cast<int>(0.080 * sr);
analysisWindowSamples = static_cast<int>(windowMs * 0.001 * sr);
// read from the history buffer size analysis +
size_t bufferSize = static_cast<size_t>(analysisWindowSamples + samples80ms + 64);
energyHistory.assign(bufferSize, 0.0f);
reset();
}
void AcousticMetrics::reset() noexcept {
std::fill(energyHistory.begin(), energyHistory.end(), 0.0f);
historyWritePos = 0;
recent50msEnergy = 0.0;
recent80msEnergy = 0.0;
totalEnergy = 0.0;
energyPeak = 0.0f;
energyPeakPos = 0;
updateCounter = 0;
d50.store(0.0f, std::memory_order_relaxed);
c50.store(0.0f, std::memory_order_relaxed);
c80.store(0.0f, std::memory_order_relaxed);
edt.store(0.0f, std::memory_order_relaxed);
}
void AcousticMetrics::processSample(float sample) noexcept {
if (energyHistory.empty()) return;
const int bufferSize = static_cast<int>(energyHistory.size());
// current sample energy ( squared )
float currentEnergy = sample * sample;
// read from the history buffer
energyHistory[historyWritePos] = currentEnergy;
// update the running sums (50 ms / 80 ms / full window)
// add the current sample, subtract the value from 50 ms ago
const int read50Pos = (historyWritePos - samples50ms + bufferSize) % bufferSize;
const int read80Pos = (historyWritePos - samples80ms + bufferSize) % bufferSize;
const int readWindowPos = (historyWritePos - analysisWindowSamples + bufferSize) % bufferSize;
recent50msEnergy += currentEnergy - energyHistory[read50Pos];
recent80msEnergy += currentEnergy - energyHistory[read80Pos];
totalEnergy += currentEnergy - energyHistory[readWindowPos];
// peak detection (EDT estimate )
if (currentEnergy > energyPeak) {
energyPeak = currentEnergy;
energyPeakPos = historyWritePos;
}
//
historyWritePos = (historyWritePos + 1) % bufferSize;
// value stable ( cumulative value 0 )
if (recent50msEnergy < 0.0) recent50msEnergy = 0.0;
if (recent80msEnergy < 0.0) recent80msEnergy = 0.0;
if (totalEnergy < 0.0) totalEnergy = 0.0;
// value interval update
if (++updateCounter >= kUpdateInterval) {
updateMetrics();
updateCounter = 0;
}
}
void AcousticMetrics::updateMetrics() noexcept {
// 50ms energy
double energy50ToInf = totalEnergy - recent50msEnergy;
if (energy50ToInf < 1e-12) energy50ToInf = 1e-12;
// 80ms energy
double energy80ToInf = totalEnergy - recent80msEnergy;
if (energy80ToInf < 1e-12) energy80ToInf = 1e-12;
// entire energy ( minimum value clipping )
double totalSafe = std::max(1e-12, totalEnergy);
// -- D50 compute (0~1 ) --
float d50val = static_cast<float>(recent50msEnergy / totalSafe);
d50val = std::min(1.0f, std::max(0.0f, d50val));
d50.store(d50val, std::memory_order_relaxed);
// -- C50 compute (dB) --
float c50val = static_cast<float>(10.0 * std::log10(recent50msEnergy / energy50ToInf));
c50val = std::min(60.0f, std::max(-60.0f, c50val));
c50.store(c50val, std::memory_order_relaxed);
// -- C80 compute (dB) --
float c80val = static_cast<float>(10.0 * std::log10(recent80msEnergy / energy80ToInf));
c80val = std::min(60.0f, std::max(-60.0f, c80val));
c80.store(c80val, std::memory_order_relaxed);
// -- EDT estimate (running) --
// after the peak, the time until energy falls to 1/10 (-10 dB decay)
// * exact EDT needs offline IR analysis; here we estimate from the peak decay time
float edtVal = 0.0f;
if (energyPeak > 1e-9f) {
// analysis peak
// scan from the peak sample until the energy reaches 1/10
const int bufferSize = static_cast<int>(energyHistory.size());
int searchStart = energyPeakPos;
float threshold = energyPeak * 0.1f; // 10dB decay
int decaySamples = 0;
for (int i = 1; i < analysisWindowSamples; ++i) {
int pos = (searchStart + i) % bufferSize;
if (energyHistory[pos] < threshold) {
decaySamples = i;
break;
}
}
edtVal = static_cast<float>(decaySamples) / static_cast<float>(sampleRate) * 6.0f;
// * 10 dB decay time x 6 ~= EDT (60 dB decay correction)
}
edt.store(edtVal, std::memory_order_relaxed);
}
// -----------------------------------------------------------------------------
// drawing: get instantaneous energy at a past time offset
// -----------------------------------------------------------------------------
// secondsAgo: how many seconds in the past to look up
// returns: the energy value at that time (squared)
//
// reads the history buffer directly for the GUI.
// out of range returns 0.
// -----------------------------------------------------------------------------
float AcousticMetrics::getEnergyAtTimeOffset(float secondsAgo) const noexcept {
if (energyHistory.empty()) return 0.0f;
const int bufferSize = static_cast<int>(energyHistory.size());
int offsetSamples = static_cast<int>(secondsAgo * static_cast<float>(sampleRate));
// clamp to range
if (offsetSamples < 0) offsetSamples = 0;
if (offsetSamples >= analysisWindowSamples) return 0.0f;
// read from the history buffer
int readPos = (historyWritePos - 1 - offsetSamples + bufferSize) % bufferSize;
return energyHistory[readPos];
}
// -----------------------------------------------------------------------------
// input activity detection: energy over the last 50 ms
// -----------------------------------------------------------------------------
// used by the GUI to hide the "measured line" when inactive.
// threshold: -60 dBFS (1e-6) energy
// -----------------------------------------------------------------------------
bool AcousticMetrics::isActive() const noexcept {
// energy over the last 50 ms determines activity
constexpr double kActivityThreshold = 1e-6; // -60dBFS
return recent50msEnergy > kActivityThreshold;
}
} // namespace FDNReverb

104
Source/DSP/AcousticMetrics.h Normal file
View file

@ -0,0 +1,104 @@
#pragma once
#include "DSPConstants.h"
#include <array>
#include <atomic>
#include <vector> // <- 1 row added
namespace FDNReverb {
// -----------------------------------------------------------------------------
// AcousticMetrics class
// -----------------------------------------------------------------------------
// Computes real-time acoustic metrics (D50, C50, C80, EDT).
//
// Principle:
// Accumulate the squared energy of the input signal in a ring buffer,
// then compare it against the energy from 50 ms / 80 ms ago,
// and compute the D50 / C50 / C80 values.
//
// Sample-rate support:
// times (ms) are converted to sample counts,
// so 44.1 kHz through 192 kHz are supported automatically.
//
// CPU:
// O(1) per-sample computation (energy accumulation)
// CPU overhead: below ~0.5%
// -----------------------------------------------------------------------------
class AcousticMetrics {
public:
AcousticMetrics() = default;
// -- initialize --
// sampleRate: sample rate (Hz)
// analysisWindowMs: analysis window (ms). 2000 ms (2 s)
void prepare(double sampleRate, float analysisWindowMs = 2000.0f);
// -- per-sample state update --
// sample: current Wet signal sample (mono)
void processSample(float sample) noexcept;
// -- value getters --
// ranges:
// D50: 0.0 ~ 1.0 ( 0.3~0.9)
// C50: -10 ~ +30 dB
// C80: -10 ~ +30 dB
// EDT: 0.0 ~ 5.0 (s)
float getD50() const noexcept { return d50.load(std::memory_order_relaxed); }
float getC50() const noexcept { return c50.load(std::memory_order_relaxed); }
float getC80() const noexcept { return c80.load(std::memory_order_relaxed); }
float getEDT() const noexcept { return edt.load(std::memory_order_relaxed); }
// --- added: expose energy history for drawing ---
// get the instantaneous energy (squared) at a time offset in the past
float getEnergyAtTimeOffset(float secondsAgo) const noexcept;
// input activity detection (energy over the last 50 ms)
bool isActive() const noexcept;
// -- reset --
void reset() noexcept;
private:
// -- compute --
void updateMetrics() noexcept;
// -- parameter --
double sampleRate{ 48000.0 };
float analysisWindowMs{ 2000.0f };
// 50ms / 80ms sample count ( sample rate depends on )
int samples50ms{ 2400 }; // @ 48kHz
int samples80ms{ 3840 }; // @ 48kHz
int analysisWindowSamples{ 96000 }; // 2000ms @ 48kHz
// -- buffer --
// energy history ( squared value )
std::vector<float> energyHistory;
int historyWritePos{ 0 };
// -- cumulative energy value --
// 50ms cumulative energy ( time )
double recent50msEnergy{ 0.0 };
// 80ms cumulative energy ( time )
double recent80msEnergy{ 0.0 };
// entire cumulative energy
double totalEnergy{ 0.0 };
// EDT : energy decay tracking
float energyPeak{ 0.0f };
int energyPeakPos{ 0 };
// -- output value (atomic for thread safety) --
std::atomic<float> d50{ 0.0f };
std::atomic<float> c50{ 0.0f };
std::atomic<float> c80{ 0.0f };
std::atomic<float> edt{ 0.0f };
// -- update --
// sample compute ,
// sample interval update
int updateCounter{ 0 };
static constexpr int kUpdateInterval = 1024; // about 21ms @ 48kHz
};
} // namespace FDNReverb

103
Source/DSP/BiquadFilters.cpp Normal file
View file

@ -0,0 +1,103 @@
#include "BiquadFilters.h"
#include "MagnitudeResponseFitter.h"
#include <JuceHeader.h>
#include <cmath>
#include <algorithm>
namespace FDNReverb {
namespace FilterDesign {
static float tanPi(float f, double fs) noexcept {
return std::tan(juce::MathConstants<float>::pi * (float)(f / fs));
}
BiquadCoeffs lowShelf(float fcHz, float gainDB, double sampleRate) {
float A = std::pow(10.f, gainDB / 40.f);
float K = tanPi(fcHz, sampleRate);
BiquadCoeffs c;
if (gainDB >= 0.f) {
float norm = 1.f / (1.f + K);
c.b0 = (1.f + A * K) * norm;
c.b1 = (A * K - 1.f) * norm;
c.b2 = 0.f;
c.a1 = (K - 1.f) * norm;
c.a2 = 0.f;
}
else {
c.b0 = (1.f + K / A) / (1.f + K);
c.b1 = (K / A - 1.f) / (1.f + K);
c.b2 = 0.f;
c.a1 = (K - 1.f) / (1.f + K);
c.a2 = 0.f;
}
return c;
}
BiquadCoeffs highShelf(float fcHz, float gainDB, double sampleRate) {
float A = std::pow(10.f, gainDB / 40.f);
float K = tanPi(fcHz, sampleRate);
BiquadCoeffs c;
if (gainDB >= 0.f) {
float norm = 1.f / (1.f + K);
c.b0 = (A + K) * norm;
c.b1 = (K - A) * norm;
c.b2 = 0.f;
c.a1 = (K - 1.f) * norm;
c.a2 = 0.f;
}
else {
float norm = 1.f / (1.f + K);
c.b0 = (1.f + A * K) * norm;
c.b1 = (A * K - 1.f) * norm;
c.b2 = 0.f;
c.a1 = (K - 1.f) * norm;
c.a2 = 0.f;
}
return c;
}
BiquadCoeffs peak(float fcHz, float gainDB, float Q, double sampleRate) {
float A = std::pow(10.f, gainDB / 40.f);
float w0 = 2.f * juce::MathConstants<float>::pi * fcHz / (float)sampleRate;
float alpha = std::sin(w0) / (2.f * Q);
float cos0 = std::cos(w0);
BiquadCoeffs c;
c.a1 = 2.f * cos0 / (1.f + alpha / A);
c.a2 = (1.f - alpha / A) / (1.f + alpha / A);
c.b0 = (1.f + alpha * A) / (1.f + alpha / A);
c.b1 = -2.f * cos0 / (1.f + alpha / A);
c.b2 = (1.f - alpha * A) / (1.f + alpha / A);
return c;
}
BiquadCoeffs highPass1st(float fcHz, double sampleRate) {
float K = tanPi(fcHz, sampleRate);
float n = 1.f + K;
BiquadCoeffs c;
c.b0 = 1.f / n; c.b1 = -1.f / n; c.b2 = 0.f;
c.a1 = (K - 1.f) / n; c.a2 = 0.f;
return c;
}
// -------------------------------------------------------------------------
// designAbsorption: MagnitudeResponseFitter
// -------------------------------------------------------------------------
// keeps the existing (UniversalEngine) helper functions,
// preserving the internal Stage-1 MRF behavior.
//
// old implementation : gain + Low/High cascade
// new implementation : Jot orthogonalizing 1 filter + LF/HF correction
// -------------------------------------------------------------------------
std::array<BiquadCoeffs, ABSO_STAGES> designAbsorption(
int delaySamples, double sampleRate,
const std::array<float, NUM_BANDS>& rt60,
float hfDamping, float lfAbsorption)
{
// MagnitudeResponseFitter processing
auto result = MagnitudeResponseFitter::design(
delaySamples, sampleRate, rt60, hfDamping, lfAbsorption);
return result.coeffs;
}
} // namespace FilterDesign
} // namespace FDNReverb

43
Source/DSP/BiquadFilters.h Normal file
View file

@ -0,0 +1,43 @@
#pragma once
#include "DSPConstants.h"
#include "../AlgorithmPresets.h"
#include <array>
namespace FDNReverb {
// -----------------------------------------------------------------------------
// Biquad helpers (Direct Form II Transposed - most robust)
// -----------------------------------------------------------------------------
struct BiquadCoeffs {
float b0{ 1.f }, b1{ 0.f }, b2{ 0.f };
float a1{ 0.f }, a2{ 0.f };
};
struct BiquadState {
float s1{ 0.f }, s2{ 0.f };
inline float tick(float x, const BiquadCoeffs& c) noexcept {
float y = c.b0 * x + s1;
s1 = c.b1 * x - c.a1 * y + s2;
s2 = c.b2 * x - c.a2 * y;
return y;
}
void reset() noexcept { s1 = s2 = 0.f; }
};
// -----------------------------------------------------------------------------
// Filter design utilities
// -----------------------------------------------------------------------------
namespace FilterDesign {
BiquadCoeffs lowShelf(float fcHz, float gainDB, double sampleRate);
BiquadCoeffs highShelf(float fcHz, float gainDB, double sampleRate);
BiquadCoeffs peak(float fcHz, float gainDB, float Q, double sampleRate);
BiquadCoeffs highPass1st(float fcHz, double sampleRate);
BiquadCoeffs allpass1st(float fcHz, double sampleRate);
// Design absorption filter cascade for delay lines
// : function internal MagnitudeResponseFitter
std::array<BiquadCoeffs, ABSO_STAGES> designAbsorption(
int delaySamples, double sampleRate,
const std::array<float, NUM_BANDS>& rt60,
float hfDamping, float lfAbsorption);
}
} // namespace FDNReverb

26
Source/DSP/DSPConstants.h Normal file
View file

@ -0,0 +1,26 @@
#pragma once
#include <array>
namespace FDNReverb {
// -- Compile-time constants ----------------------------------------------------
static constexpr int FDN_N = 8; // FDN order (channels; legacy definition kept for reference)
static constexpr int SAPF_STAGES = 3; // allpass stages per delay line
static constexpr int ABSO_STAGES = 3; // Stage 1: Jot first-order + LF/HF correction
static constexpr int ER_TAPS = 16; // early-reflection FIR taps
// Stage 2 (Valimaki-Liski cumulative GEQ) stages:
// 10: 10-band GEQ (interaction matrix + WLS)
//
// important design notes:
// - the mid-band gain (midGain) of GEQ band 0 is absorbed into the b0/b1/b2 coefficients;
// no separate DC gain stage is needed to avoid DC coloration,
// just a single gain.
// - LF Absorption / HF Damping are applied directly as GEQ target dB,
// fully independent of each other.
// - targets are clamped to 0 dB or below, mathematically guaranteeing loop gain <= 1.
static constexpr int ABSO_STAGES_S2 = 10;
// Mutually-prime base delays (samples @ 48 kHz), log-distributed 30-130 ms
static constexpr std::array<int, FDN_N> BASE_PRIMES_48K = {
1451, 1693, 1979, 2311, 2683, 3067, 3491, 3923
};
} // namespace FDNReverb

134
Source/DSP/DelayMemory.h Normal file
View file

@ -0,0 +1,134 @@
#pragma once
#include <vector>
#include <cmath>
#include <algorithm>
#include <cstdint>
namespace FDNReverb {
// -----------------------------------------------------------------------------
// memory pool (Single-Large Buffer)
// -----------------------------------------------------------------------------
class DelayMemoryPool {
public:
void allocate(size_t totalSamples) {
buffer.assign(totalSamples, 0.0f);
allocOffset = 0;
}
// pointer sized up to the next power of two (also outputs an index mask)
float* requestMemory(size_t samplesNeeded, int& outMask) {
size_t powerOfTwoSize = 1;
while (powerOfTwoSize < samplesNeeded) powerOfTwoSize *= 2;
if (allocOffset + powerOfTwoSize > buffer.size()) return nullptr;
float* ptr = buffer.data() + allocOffset;
outMask = static_cast<int>(powerOfTwoSize - 1);
allocOffset += powerOfTwoSize;
return ptr;
}
void clear() { std::fill(buffer.begin(), buffer.end(), 0.0f); }
private:
std::vector<float> buffer;
size_t allocOffset{ 0 };
};
// -----------------------------------------------------------------------------
// interpolation
// -----------------------------------------------------------------------------
class LinearDelayLine {
public:
void init(float* memory, int bitmask) {
buffer = memory;
mask = bitmask;
writeIndex = 0;
}
// linear interpolation ( high band natural Air Absorption )
inline float read(float delayInSamples) const noexcept {
int id = static_cast<int>(delayInSamples);
float frac = delayInSamples - static_cast<float>(id);
// bitwise ops undefined behavior completely , uint32_t
uint32_t uWrite = static_cast<uint32_t>(writeIndex);
uint32_t uId = static_cast<uint32_t>(id);
uint32_t uMask = static_cast<uint32_t>(mask);
int readIdx1 = static_cast<int>((uWrite - uId) & uMask);
int readIdx2 = static_cast<int>((uWrite - uId - 1) & uMask);
return buffer[readIdx1] + frac * (buffer[readIdx2] - buffer[readIdx1]);
}
inline void write(float input) noexcept {
buffer[writeIndex] = input;
writeIndex = (writeIndex + 1) & mask;
}
private:
float* buffer{ nullptr };
int mask{ 0 };
int writeIndex{ 0 };
};
// -----------------------------------------------------------------------------
// Thiran allpass interpolation (preserves the phase response)
// linear interpolation would dull high-band decay (sinc(pi*f) rolloff), so use a Thiran allpass
// which keeps |H(w)| = 1, preserving high-band clarity in the FDN feedback loops.
// -----------------------------------------------------------------------------
class ThiranDelayLine {
public:
void init(float* memory, int bitmask) {
buffer = memory;
mask = bitmask;
writeIndex = 0;
thiranX1 = 0.0f;
thiranY1 = 0.0f;
}
void resetState() noexcept {
thiranX1 = 0.0f;
thiranY1 = 0.0f;
}
// Thiran first-order allpass: y[n] = a*x[n] + x[n-1] - a*y[n-1]
// a = (1-D)/(1+D), D = fractional delay
inline float read(float delayInSamples) noexcept {
int id = static_cast<int>(delayInSamples);
float frac = delayInSamples - static_cast<float>(id);
// clamp below to avoid instability as frac->0, a->1
frac = std::max(frac, 0.1f);
const float a = (1.0f - frac) / (1.0f + frac);
uint32_t uWrite = static_cast<uint32_t>(writeIndex);
uint32_t uId = static_cast<uint32_t>(id);
uint32_t uMask = static_cast<uint32_t>(mask);
float xn = buffer[static_cast<int>((uWrite - uId) & uMask)];
float yn = a * xn + thiranX1 - a * thiranY1;
thiranX1 = xn;
thiranY1 = yn;
return yn;
}
inline void write(float input) noexcept {
buffer[writeIndex] = input;
writeIndex = (writeIndex + 1) & mask;
}
private:
float* buffer{ nullptr };
int mask{ 0 };
int writeIndex{ 0 };
float thiranX1{ 0.0f };
float thiranY1{ 0.0f };
};
} // namespace FDNReverb

66
Source/DSP/EarlyReflections.cpp Normal file
View file

@ -0,0 +1,66 @@
#include "EarlyReflections.h"
namespace FDNReverb {
void EarlyReflections::prepare(const juce::dsp::ProcessSpec& spec) {
int maxSamples = (int)(0.7 * spec.sampleRate) + 8;
juce::dsp::ProcessSpec mono = spec;
mono.numChannels = 1;
buf.prepare(mono);
buf.setMaximumDelayInSamples(maxSamples);
erHPCoeffs = FilterDesign::highPass1st(80.f, spec.sampleRate);
float K = std::tan(juce::MathConstants<float>::pi * 6000.f / (float)spec.sampleRate);
erLPCoeffs.b0 = K / (1.f + K);
erLPCoeffs.b1 = erLPCoeffs.b0;
erLPCoeffs.b2 = 0.f;
erLPCoeffs.a1 = (K - 1.f) / (K + 1.f);
erLPCoeffs.a2 = 0.f;
}
void EarlyReflections::buildTaps(const AlgorithmPreset& preset, float roomSizeScale, double sampleRate) {
float erEnergy50 = preset.acoustics.d50[4];
float V = preset.volumeM3 > 0.f ? preset.volumeM3 : 10.f;
float mixTimeMs = std::min(0.0117f * V + 50.1f, 150.f);
float span = mixTimeMs * roomSizeScale;
for (int i = 0; i < ER_TAPS; ++i) {
float t01 = static_cast<float>(i + 1) / static_cast<float>(ER_TAPS);
float delMs = span * std::pow(t01, 1.5f);
taps[i].delaySamples = delMs * 0.001f * (float)sampleRate;
float rt60m = preset.acoustics.rt60[4];
float amp = std::exp(-6.9f * delMs * 0.001f / rt60m);
float factor = (i < ER_TAPS / 2) ? std::sqrt(erEnergy50) : std::sqrt(1.f - erEnergy50);
amp *= factor * std::sqrt(2.f / ER_TAPS);
float pan = (i % 3 == 0) ? -0.707f : ((i % 3 == 1) ? 0.707f : 0.0f);
taps[i].gainL = amp * std::sqrt(0.5f - 0.5f * pan);
taps[i].gainR = amp * std::sqrt(0.5f + 0.5f * pan);
}
}
void EarlyReflections::setPreDelay(float ms, double sampleRate) noexcept {
preDelaySamples = juce::roundToInt(ms * 0.001 * sampleRate);
}
std::pair<float, float> EarlyReflections::tick(float mono) noexcept {
buf.pushSample(0, mono);
float L = 0.f, R = 0.f;
for (const auto& t : taps) {
float d = buf.popSample(0, t.delaySamples + preDelaySamples, false);
L += t.gainL * d;
R += t.gainR * d;
}
L = erHPL.tick(L, erHPCoeffs);
R = erHPR.tick(R, erHPCoeffs);
return { L, R };
}
void EarlyReflections::reset() noexcept {
buf.reset();
erHPL.reset(); erHPR.reset();
erLPL.reset(); erLPR.reset();
}
} // namespace FDNReverb

32
Source/DSP/EarlyReflections.h Normal file
View file

@ -0,0 +1,32 @@
#pragma once
#include <JuceHeader.h>
#include "DSPConstants.h"
#include "BiquadFilters.h"
namespace FDNReverb {
struct ERTap {
float delaySamples{ 0.f };
float gainL{ 0.f };
float gainR{ 0.f };
};
class EarlyReflections {
public:
void prepare(const juce::dsp::ProcessSpec& spec);
void buildTaps(const AlgorithmPreset& preset, float roomSizeScale, double sampleRate);
void setPreDelay(float ms, double sampleRate) noexcept;
std::pair<float, float> tick(float mono) noexcept;
void reset() noexcept;
private:
juce::dsp::DelayLine<float, juce::dsp::DelayLineInterpolationTypes::Lagrange3rd> buf;
std::array<ERTap, ER_TAPS> taps;
int preDelaySamples{ 0 };
BiquadCoeffs erHPCoeffs, erLPCoeffs;
BiquadState erHPL, erHPR, erLPL, erLPR;
};
} // namespace FDNReverb

388
Source/DSP/MagnitudeResponseFitter.cpp Normal file
View file

@ -0,0 +1,388 @@
#include "MagnitudeResponseFitter.h"
#include <JuceHeader.h>
#include <cmath>
#include <algorithm>
#include <complex>
namespace FDNReverb {
// -----------------------------------------------------------------------------
// static
// -----------------------------------------------------------------------------
std::array<std::array<double, NUM_BANDS>, NUM_BANDS> MagnitudeResponseFitter::cachedB;
std::array<std::array<double, NUM_BANDS>, NUM_BANDS> MagnitudeResponseFitter::cachedBtWB;
std::array<double, NUM_BANDS> MagnitudeResponseFitter::cachedW;
double MagnitudeResponseFitter::cachedSampleRate = 0.0;
bool MagnitudeResponseFitter::cacheValid = false;
// -----------------------------------------------------------------------------
// band Q value ( band : Q ~ sqrt2 / (2^(1/2) - 2^(-1/2)) ~ 1.414)
// -----------------------------------------------------------------------------
static const std::array<float, NUM_BANDS> kBandQs = {
1.7f, // 31.25 Hz (: Q rise )
1.414f, // 62.5 Hz
1.414f, // 125 Hz
1.414f, // 250 Hz
1.414f, // 500 Hz
1.414f, // 1 kHz
1.414f, // 2 kHz
1.414f, // 4 kHz
1.414f, // 8 kHz
1.7f // 16 kHz (: Q rise )
};
const std::array<float, NUM_BANDS>& MagnitudeResponseFitter::getBandQs() noexcept {
return kBandQs;
}
// -----------------------------------------------------------------------------
// Stage 1 ( existing )
// -----------------------------------------------------------------------------
float MagnitudeResponseFitter::t60ToLoopGain(float t60Seconds, int delaySamples, double sampleRate) noexcept {
float t60Safe = std::max(0.01f, t60Seconds);
float exponent = -3.0f * static_cast<float>(delaySamples) / (static_cast<float>(sampleRate) * t60Safe);
return std::pow(10.0f, exponent);
}
float MagnitudeResponseFitter::computeJotPole(float gDC, float alphaRatio) noexcept {
float alphaSafe = juce::jlimit(0.05f, 20.0f, alphaRatio);
float gDCSafe = juce::jlimit(1e-6f, 0.99999f, gDC);
constexpr float kLn10Over4 = 0.5756462732485f;
float log10g = std::log10(gDCSafe);
float alphaSqInv = 1.0f / (alphaSafe * alphaSafe);
float pole = kLn10Over4 * log10g * (1.0f - alphaSqInv);
return juce::jlimit(-0.98f, 0.98f, pole);
}
BiquadCoeffs MagnitudeResponseFitter::orthogonalizedFirstOrderToBiquad(float gain, float pole) noexcept {
BiquadCoeffs c;
c.b0 = gain * (1.0f - pole);
c.b1 = 0.0f;
c.b2 = 0.0f;
c.a1 = -pole;
c.a2 = 0.0f;
return c;
}
float MagnitudeResponseFitter::getT60AtDC(const std::array<float, NUM_BANDS>& rt60) noexcept {
return (rt60[0] + rt60[1]) * 0.5f;
}
float MagnitudeResponseFitter::getT60AtNyquist(const std::array<float, NUM_BANDS>& rt60, double sampleRate) noexcept {
if (sampleRate <= 50000.0) {
return rt60[9];
}
else {
return (rt60[8] + rt60[9]) * 0.5f;
}
}
// -----------------------------------------------------------------------------
// Stage 1 main design function ( existing )
// -----------------------------------------------------------------------------
MagnitudeResponseFitter::DesignResult MagnitudeResponseFitter::design(
int delaySamples,
double sampleRate,
const std::array<float, NUM_BANDS>& rt60,
float hfDamping,
float lfAbsorption)
{
DesignResult result;
float t60DC = std::max(0.01f, getT60AtDC(rt60));
float t60Nyq = std::max(0.01f, getT60AtNyquist(rt60, sampleRate));
float gDC = t60ToLoopGain(t60DC, delaySamples, sampleRate);
float gNyq = t60ToLoopGain(t60Nyq, delaySamples, sampleRate);
float alpha = t60Nyq / t60DC;
float pole = computeJotPole(gDC, alpha);
result.coeffs[0] = orthogonalizedFirstOrderToBiquad(gDC, pole);
float lfShelfDB = -lfAbsorption * 3.0f;
result.coeffs[1] = FilterDesign::lowShelf(150.0f, lfShelfDB, sampleRate);
float hfShelfDB = -hfDamping * 6.0f;
result.coeffs[2] = FilterDesign::highShelf(4000.0f, hfShelfDB, sampleRate);
result.dcGain = gDC;
result.nyquistGain = gNyq;
result.pole = pole;
return result;
}
// -----------------------------------------------------------------------------
// Stage 2 : Biquad peak filter
// -----------------------------------------------------------------------------
BiquadCoeffs MagnitudeResponseFitter::designSymmetricPeakBiquad(
float fcHz, float gainDB, float Q, double sampleRate) noexcept
{
float fcSafe = juce::jlimit(10.0f, static_cast<float>(sampleRate) * 0.49f, fcHz);
float A = std::pow(10.0f, gainDB / 40.0f);
float w0 = 2.0f * juce::MathConstants<float>::pi * fcSafe / static_cast<float>(sampleRate);
float cosW0 = std::cos(w0);
float sinW0 = std::sin(w0);
float alpha = sinW0 / (2.0f * std::max(0.1f, Q));
float a0 = 1.0f + alpha / A;
BiquadCoeffs c;
c.b0 = (1.0f + alpha * A) / a0;
c.b1 = -2.0f * cosW0 / a0;
c.b2 = (1.0f - alpha * A) / a0;
c.a1 = -2.0f * cosW0 / a0;
c.a2 = (1.0f - alpha / A) / a0;
return c;
}
// -----------------------------------------------------------------------------
// Stage 2 : Biquad magnitude response (dB) compute
// -----------------------------------------------------------------------------
float MagnitudeResponseFitter::biquadMagnitudeDB(
const BiquadCoeffs& c, float fEval, double sampleRate) noexcept
{
double w = 2.0 * juce::MathConstants<double>::pi * fEval / sampleRate;
double cosW = std::cos(w);
double sinW = std::sin(w);
double cos2W = std::cos(2.0 * w);
double sin2W = std::sin(2.0 * w);
double bRe = c.b0 + c.b1 * cosW + c.b2 * cos2W;
double bIm = -c.b1 * sinW - c.b2 * sin2W;
double aRe = 1.0 + c.a1 * cosW + c.a2 * cos2W;
double aIm = -c.a1 * sinW - c.a2 * sin2W;
double bMag2 = bRe * bRe + bIm * bIm;
double aMag2 = aRe * aRe + aIm * aIm;
double mag2 = bMag2 / std::max(1e-30, aMag2);
return static_cast<float>(10.0 * std::log10(std::max(1e-30, mag2)));
}
// -----------------------------------------------------------------------------
// Stage 2 : 10x10 LDLT decomposition solver
// -----------------------------------------------------------------------------
void MagnitudeResponseFitter::solveLDLT10(
const std::array<std::array<double, NUM_BANDS>, NUM_BANDS>& A,
const std::array<double, NUM_BANDS>& b,
std::array<double, NUM_BANDS>& x) noexcept
{
constexpr int N = NUM_BANDS;
double L[N][N] = { 0 };
double D[N] = { 0 };
for (int i = 0; i < N; ++i) L[i][i] = 1.0;
for (int j = 0; j < N; ++j) {
double sum = A[j][j];
for (int k = 0; k < j; ++k) {
sum -= L[j][k] * L[j][k] * D[k];
}
D[j] = sum;
if (std::abs(D[j]) < 1e-12) {
D[j] = (D[j] < 0.0 ? -1e-12 : 1e-12);
}
for (int i = j + 1; i < N; ++i) {
double s = A[i][j];
for (int k = 0; k < j; ++k) {
s -= L[i][k] * L[j][k] * D[k];
}
L[i][j] = s / D[j];
}
}
double z[N];
for (int i = 0; i < N; ++i) {
double s = b[i];
for (int k = 0; k < i; ++k) s -= L[i][k] * z[k];
z[i] = s;
}
double y[N];
for (int i = 0; i < N; ++i) y[i] = z[i] / D[i];
for (int i = N - 1; i >= 0; --i) {
double s = y[i];
for (int k = i + 1; k < N; ++k) s -= L[k][i] * x[k];
x[i] = s;
}
}
// -----------------------------------------------------------------------------
// Stage 2 : Biquad coefficient linear gain absorption
// -----------------------------------------------------------------------------
// H(z) = (b0 + b1.z^{-1} + b2.z^{-2}) / (1 + a1.z^{-1} + a2.z^{-2})
//
// frequency amplitude linearGain , (b0, b1, b2) linearGain
// . mathematically independent DC color apply completely .
BiquadCoeffs MagnitudeResponseFitter::absorbGainIntoBiquad(
const BiquadCoeffs& c, float linearGain) noexcept
{
BiquadCoeffs result = c;
result.b0 *= linearGain;
result.b1 *= linearGain;
result.b2 *= linearGain;
return result;
}
// -----------------------------------------------------------------------------
// Stage 2: Interaction Matrix before compute
// -----------------------------------------------------------------------------
void MagnitudeResponseFitter::precomputeInteractionMatrix(double sampleRate) {
if (cacheValid && std::abs(cachedSampleRate - sampleRate) < 0.5) {
return;
}
constexpr int N = NUM_BANDS;
constexpr float kProbeGainDB = 1.0f;
for (int j = 0; j < N; ++j) {
BiquadCoeffs c = designSymmetricPeakBiquad(
BAND_FREQ[j], kProbeGainDB, kBandQs[j], sampleRate);
for (int i = 0; i < N; ++i) {
float dB = biquadMagnitudeDB(c, BAND_FREQ[i], sampleRate);
cachedB[i][j] = static_cast<double>(dB);
}
}
const std::array<double, NUM_BANDS> weights = {
0.5, // 31.25 Hz
0.7, // 62.5 Hz
0.85, // 125 Hz
1.0, // 250 Hz
1.0, // 500 Hz
1.0, // 1 kHz
1.0, // 2 kHz
1.0, // 4 kHz
0.85, // 8 kHz
0.6 // 16 kHz
};
for (int i = 0; i < N; ++i) cachedW[i] = weights[i];
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
double s = 0.0;
for (int k = 0; k < N; ++k) {
s += cachedB[k][i] * cachedW[k] * cachedB[k][j];
}
cachedBtWB[i][j] = s;
}
}
constexpr double kRidge = 1e-4;
for (int i = 0; i < N; ++i) cachedBtWB[i][i] += kRidge;
cachedSampleRate = sampleRate;
cacheValid = true;
}
// -----------------------------------------------------------------------------
// Stage 2c: main design function ( fix )
// -----------------------------------------------------------------------------
// :
// 1. band target dB compute (T60 dB )
// t[i] = -60 . m / (fs . T60[i])
// 2. LF/HF correction target dB directly
// 3. target dB 0 below clamp -> loop gain <= 1 guarantee
// 4. mid-band gain midGain (band 4 = 500Hz)
// midGain = 10^(midDb/20)
// 5. dB WLS
// g_cmd = (B^T.W.B)^(-1).B^T.W.t_residual
// 6. g_cmd[j] dB Biquad coefficient
// 7. band 0 coefficient midGain absorption
// -> independent DC color apply not needed
MagnitudeResponseFitter::DesignResultStage2 MagnitudeResponseFitter::designStage2(
int delaySamples,
double sampleRate,
const std::array<float, NUM_BANDS>& rt60,
float hfDamping,
float lfAbsorption)
{
precomputeInteractionMatrix(sampleRate);
DesignResultStage2 result;
constexpr int N = NUM_BANDS;
const float fs = static_cast<float>(sampleRate);
const float m = static_cast<float>(delaySamples);
// -- Step 1: band loop 1 gain dB target --
std::array<float, NUM_BANDS> targetDb;
for (int i = 0; i < N; ++i) {
float t60Safe = std::max(0.01f, rt60[i]);
targetDb[i] = -60.0f * m / (fs * t60Safe);
}
// -- Step 2: LF/HF correction target dB --
// LF Absorption: low band (31Hz, 62Hz, 125Hz) added decay
// lfAbsorption=0 -> correction , =1 -> -3dB added decay
targetDb[0] += -lfAbsorption * 3.0f;
targetDb[1] += -lfAbsorption * 2.5f;
targetDb[2] += -lfAbsorption * 1.5f;
// HF Damping: high band (4kHz, 8kHz, 16kHz) added decay
// hfDamping=0 -> correction , =1 -> -6dB added decay
targetDb[7] += -hfDamping * 3.0f;
targetDb[8] += -hfDamping * 5.0f;
targetDb[9] += -hfDamping * 6.0f;
// -- Step 3: target dB 0 below clamp --
// loop gain <= 1 mathematically guarantee safe
for (int i = 0; i < N; ++i) {
targetDb[i] = std::min(targetDb[i], 0.0f);
// decay precision influence below (-60dB/loop)
targetDb[i] = std::max(targetDb[i], -60.0f);
result.targetDb[i] = targetDb[i];
}
// -- Step 4: mid-band gain midGain (band 4 = 500Hz) --
float midDb = targetDb[4];
float midGainLinear = std::pow(10.0f, midDb / 20.0f);
result.midGainAbsorbed = midGainLinear;
// dB: mid-band deviation (GEQ frequency response )
std::array<double, NUM_BANDS> residualDb;
for (int i = 0; i < N; ++i) {
residualDb[i] = static_cast<double>(targetDb[i] - midDb);
}
// -- Step 5: WLS GEQ coefficient --
std::array<double, NUM_BANDS> rhs;
for (int j = 0; j < N; ++j) {
double s = 0.0;
for (int k = 0; k < N; ++k) {
s += cachedB[k][j] * cachedW[k] * residualDb[k];
}
rhs[j] = s;
}
std::array<double, NUM_BANDS> gCmd;
solveLDLT10(cachedBtWB, rhs, gCmd);
// -- Step 6: g_cmd[j] dB Biquad coefficient --
// safe range clamp (+/-18 dB )
for (int j = 0; j < N; ++j) {
float gDb = static_cast<float>(juce::jlimit(-18.0, 18.0, gCmd[j]));
result.commandDb[j] = gDb;
result.geqStages[j] = designSymmetricPeakBiquad(
BAND_FREQ[j], gDb, kBandQs[j], sampleRate);
}
// -- Step 7: band 0 coefficient midGain absorption --
// independent DC color apply not needed ,
// filter cascade entire loop gain exact WLS .
result.geqStages[0] = absorbGainIntoBiquad(result.geqStages[0], midGainLinear);
return result;
}
} // namespace FDNReverb

130
Source/DSP/MagnitudeResponseFitter.h Normal file
View file

@ -0,0 +1,130 @@
#pragma once
#include "DSPConstants.h"
#include "BiquadFilters.h"
#include "../AlgorithmPresets.h"
#include <array>
namespace FDNReverb {
// -----------------------------------------------------------------------------
// MagnitudeResponseFitter
// -----------------------------------------------------------------------------
// designs the 10-band RT60 absorption filters for the FDN.
//
// design modes :
// Stage 1 (Jot first-order orthogonalizing):
// Jot-Chaigne (AES Preprint 3030, 1991) first-order orthogonalizing filters.
// matched at DC and Nyquist with 2 design points.
//
// Stage 2c (Valimaki-Liski cumulative GEQ):
// Valimaki & Liski (IEEE SPL 2017) Interaction Matrix + WLS
// exact fit across the 10 bands.
//
// safety guarantee :
// - targets are clamped to 0 dB or below -> loop gain <= 1 is guaranteed
// - band 0 midGain and b0/b1/b2 are absorbed into the applied filter
// - LF/HF corrections are independent GEQ targets in dB
//
// important :
// - per-band decay in dB: -60*m / (fs*T60)
// avoids the "2 kHz T60 assumption" of Schlecht-Habets (DAFx-17)
// - design runs offline (message thread); the resulting Biquad coefficients
// are used on the audio thread
// -----------------------------------------------------------------------------
class MagnitudeResponseFitter {
public:
enum class DesignMode {
Stage1_Jot1stOrder, // Jot first-order orthogonalizing (2 pts: DC/Nyquist)
Stage2_BiquadGEQ // Valimaki-Liski cumulative GEQ (exact at 10 bands)
};
// -------------------------------------------------------------------------
// Stage 1 design result (existing)
// -------------------------------------------------------------------------
// ABSO_STAGES = 3 Biquads:
// coeffs[0] = gain (Jot first-order orthogonalizing filter, Biquad form)
// coeffs[1] = low-band correction (Low Shelf, LF Absorption)
// coeffs[2] = high-band correction (High Shelf, HF Damping)
struct DesignResult {
std::array<BiquadCoeffs, ABSO_STAGES> coeffs;
float dcGain{ 1.0f };
float nyquistGain{ 1.0f };
float pole{ 0.0f };
};
// -------------------------------------------------------------------------
// Stage 2c design result
// -------------------------------------------------------------------------
// 10-band GEQ:
// geqStages[0] = band 0 (31.25 Hz), midGain absorbed into the coefficient
// geqStages[1..9] = bands 1-9 (62.5 Hz - 16 kHz), GEQ
//
// filter chain: geqStages[0] -> geqStages[1] -> ... -> geqStages[9]
// no separate midGain stage is needed (absorbed into band 0).
struct DesignResultStage2 {
std::array<BiquadCoeffs, NUM_BANDS> geqStages; // 10-band GEQ
// visualization
std::array<float, NUM_BANDS> targetDb; // per-band target dB (after clamping)
std::array<float, NUM_BANDS> commandDb; // WLS-solved command dB
float midGainAbsorbed{ 1.0f }; // midGain absorbed into band 0
};
// -------------------------------------------------------------------------
// Stage 1 design function (existing)
// -------------------------------------------------------------------------
static DesignResult design(
int delaySamples,
double sampleRate,
const std::array<float, NUM_BANDS>& rt60,
float hfDamping,
float lfAbsorption);
// -------------------------------------------------------------------------
// Stage 2c design function
// -------------------------------------------------------------------------
static DesignResultStage2 designStage2(
int delaySamples,
double sampleRate,
const std::array<float, NUM_BANDS>& rt60,
float hfDamping,
float lfAbsorption);
// -------------------------------------------------------------------------
// precompute the interaction matrix once (per sample rate)
// -------------------------------------------------------------------------
static void precomputeInteractionMatrix(double sampleRate);
static double getCachedSampleRate() noexcept { return cachedSampleRate; }
private:
// -- Stage 1 --
static float t60ToLoopGain(float t60Seconds, int delaySamples, double sampleRate) noexcept;
static float computeJotPole(float gDC, float alphaRatio) noexcept;
static BiquadCoeffs orthogonalizedFirstOrderToBiquad(float gain, float pole) noexcept;
static float getT60AtDC(const std::array<float, NUM_BANDS>& rt60) noexcept;
static float getT60AtNyquist(const std::array<float, NUM_BANDS>& rt60, double sampleRate) noexcept;
// -- Stage 2 --
static BiquadCoeffs designSymmetricPeakBiquad(
float fcHz, float gainDB, float Q, double sampleRate) noexcept;
static const std::array<float, NUM_BANDS>& getBandFreqs() noexcept { return BAND_FREQ; }
static const std::array<float, NUM_BANDS>& getBandQs() noexcept;
static float biquadMagnitudeDB(const BiquadCoeffs& c, float fEval, double sampleRate) noexcept;
static void solveLDLT10(
const std::array<std::array<double, NUM_BANDS>, NUM_BANDS>& A,
const std::array<double, NUM_BANDS>& b,
std::array<double, NUM_BANDS>& x) noexcept;
// absorb the entire DC gain of the Biquad (b0, b1, b2) into a gain
// so an independent DC gain can be applied to the filter mathematically
static BiquadCoeffs absorbGainIntoBiquad(const BiquadCoeffs& c, float linearGain) noexcept;
// -- Stage 2 static --
static std::array<std::array<double, NUM_BANDS>, NUM_BANDS> cachedB;
static std::array<std::array<double, NUM_BANDS>, NUM_BANDS> cachedBtWB;
static std::array<double, NUM_BANDS> cachedW;
static double cachedSampleRate;
static bool cacheValid;
};
} // namespace FDNReverb

148
Source/DSP/OutputEQ.h Normal file
View file

@ -0,0 +1,148 @@
#pragma once
#include <cmath>
#include <algorithm>
namespace FDNReverb {
// -----------------------------------------------------------------------------
// OutputEQ: Wet output stage Lo/Hi Cut (Linkwitz-Riley 12dB/oct)
// -----------------------------------------------------------------------------
// design rationale:
// - 1 IIR (6dB/oct) x 2 cascade = 12dB/oct
// - Linkwitz-Riley topology: 2nd-order phase alignment
// - keeps the reverb sounding musical
//
// filter equation (1 IIR):
// HPF: y[n] = R . (y[n-1] + x[n] - x[n-1])
// LPF: y[n] = (1 - R) . x[n] + R . y[n-1]
// where R = exp(-2pi.fc/fs)
//
// real-time safety :
// - no allocation at all
// - per-sample cost: HPF 8 ops + LPF 6 ops (L/R combined)
// - coefficients updated per block (no zipper noise, no SmoothedValue needed)
//
// bypass :
// - Lo Cut below 20 Hz -> HPF fully bypassed
// - Hi Cut above 20 kHz -> LPF fully bypassed
// both bypasses are per-block coefficient updates, so CPU use is trivial.
// -----------------------------------------------------------------------------
class OutputEQ {
public:
OutputEQ() = default;
void prepare(double sampleRate) noexcept {
fs = sampleRate;
reset();
setLoCutHz(20.0f);
setHiCutHz(20000.0f);
}
void reset() noexcept {
// HPF state (two stages per channel, L/R)
hpfX1_L_1 = hpfY1_L_1 = 0.0f;
hpfX1_L_2 = hpfY1_L_2 = 0.0f;
hpfX1_R_1 = hpfY1_R_1 = 0.0f;
hpfX1_R_2 = hpfY1_R_2 = 0.0f;
// LPF state (two stages per channel, L/R)
lpfY1_L_1 = 0.0f;
lpfY1_L_2 = 0.0f;
lpfY1_R_1 = 0.0f;
lpfY1_R_2 = 0.0f;
}
// --- parameter setters (called per block) ---
void setLoCutHz(float fcHz) noexcept {
currentLoCutHz = fcHz;
// bypass below 20 Hz (skip R computation)
if (fcHz <= 20.0f) {
loCutActive = false;
return;
}
loCutActive = true;
constexpr float twoPi = 6.28318530718f;
const float clamped = std::clamp(fcHz, 20.0f, 500.0f);
loCutR = std::exp(-twoPi * clamped / static_cast<float>(fs));
}
void setHiCutHz(float fcHz) noexcept {
currentHiCutHz = fcHz;
// bypass above 20 kHz
const float nyquist = static_cast<float>(fs) * 0.45f;
const float clamped = std::clamp(fcHz, 1000.0f, std::min(20000.0f, nyquist));
if (fcHz >= 20000.0f) {
hiCutActive = false;
return;
}
hiCutActive = true;
constexpr float twoPi = 6.28318530718f;
hiCutR = std::exp(-twoPi * clamped / static_cast<float>(fs));
}
// --- per-sample processing (L/R interleaved) ---
inline void process(float& l, float& r) noexcept {
// -- Lo Cut: 1 HPF x 2 cascade --
if (loCutActive) {
// L stage 1
const float l_in = l;
const float l_1 = loCutR * (hpfY1_L_1 + l_in - hpfX1_L_1);
hpfX1_L_1 = l_in;
hpfY1_L_1 = l_1;
// L stage 2
const float l_2 = loCutR * (hpfY1_L_2 + l_1 - hpfX1_L_2);
hpfX1_L_2 = l_1;
hpfY1_L_2 = l_2;
l = l_2;
// R stage 1
const float r_in = r;
const float r_1 = loCutR * (hpfY1_R_1 + r_in - hpfX1_R_1);
hpfX1_R_1 = r_in;
hpfY1_R_1 = r_1;
// R stage 2
const float r_2 = loCutR * (hpfY1_R_2 + r_1 - hpfX1_R_2);
hpfX1_R_2 = r_1;
hpfY1_R_2 = r_2;
r = r_2;
}
// -- Hi Cut: 1 LPF x 2 cascade --
if (hiCutActive) {
const float oneMinusR = 1.0f - hiCutR;
// L stage 1
lpfY1_L_1 = oneMinusR * l + hiCutR * lpfY1_L_1;
// L stage 2
lpfY1_L_2 = oneMinusR * lpfY1_L_1 + hiCutR * lpfY1_L_2;
l = lpfY1_L_2;
// R stage 1
lpfY1_R_1 = oneMinusR * r + hiCutR * lpfY1_R_1;
// R stage 2
lpfY1_R_2 = oneMinusR * lpfY1_R_1 + hiCutR * lpfY1_R_2;
r = lpfY1_R_2;
}
}
float getCurrentLoCutHz() const noexcept { return currentLoCutHz; }
float getCurrentHiCutHz() const noexcept { return currentHiCutHz; }
private:
double fs{ 48000.0 };
// -- Lo Cut (HPF) --
bool loCutActive{ false };
float loCutR{ 0.0f };
float currentLoCutHz{ 20.0f };
float hpfX1_L_1{}, hpfY1_L_1{}, hpfX1_L_2{}, hpfY1_L_2{};
float hpfX1_R_1{}, hpfY1_R_1{}, hpfX1_R_2{}, hpfY1_R_2{};
// -- Hi Cut (LPF) --
bool hiCutActive{ false };
float hiCutR{ 0.0f };
float currentHiCutHz{ 20000.0f };
float lpfY1_L_1{}, lpfY1_L_2{};
float lpfY1_R_1{}, lpfY1_R_2{};
};
} // namespace FDNReverb

79
Source/DSP/OutputLimiter.h Normal file
View file

@ -0,0 +1,79 @@
#pragma once
#include <cmath>
#include <algorithm>
namespace FDNReverb {
// -----------------------------------------------------------------------------
// OutputLimiter: safe output stage (true-peak limiter)
// -----------------------------------------------------------------------------
// design rationale :
// - parameter values are chosen conservatively for safety
// - Threshold = -0.5 dBFS (~0.944): suppresses peaks before the DAW limiter
// - Look-ahead: introduces plugin latency
// - Attack: 0.5 ms (peak-based)
// - Release: 50 ms (prevents unnatural pumping)
//
// real-time safety :
// - allocation: once in prepare(), never in processBlock
// - per-sample gain: one comparison against targetGain, SIMD-friendly
// - floating-point math: no branches or transcendental functions
//
// - layout: output stage of UniversalEngine::processBlock()
// (after Dry/Wet mix, before the stereo output)
// -----------------------------------------------------------------------------
class OutputLimiter {
public:
OutputLimiter() = default;
// --- sample-rate dependent coefficient computation ---
void prepare(double sampleRate) noexcept {
fs = sampleRate;
// 1 path filter coefficient : y[n] = y[n-1] + coeff * (x[n] - y[n-1])
// coeff = 1 - exp(-T / tau) where T = 1/fs, tau = time constant
attackCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * 0.0005f)); // 0.5ms
releaseCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * 0.050f)); // 50ms
reset();
}
void reset() noexcept {
currentGain = 1.0f;
}
// --- per-sample processing (called from the audio thread) ---
inline void process(float& l, float& r) noexcept {
// peak detection (max of L/R levels)
const float absL = std::abs(l);
const float absR = std::abs(r);
const float peak = std::max(absL, absR);
// Threshold: -0.5 dBFS ~ 0.944
// compute target gain from the signal
constexpr float threshold = 0.944f;
// target gain :
// peak <= threshold -> 1.0 (no reduction needed)
// peak > threshold -> threshold/peak (pull signal to threshold)
const float targetGain = (peak > threshold) ? (threshold / peak) : 1.0f;
// Attack/Release envelope
// when targetGain < currentGain (gain must decrease): attack
// when targetGain > currentGain (gain recovers): release
// so peaks are suppressed smoothly
const float coeff = (targetGain < currentGain) ? attackCoeff : releaseCoeff;
currentGain += (targetGain - currentGain) * coeff;
// apply the same gain to L/R to preserve the stereo image
l *= currentGain;
r *= currentGain;
}
private:
double fs{ 44100.0 };
float attackCoeff{ 0.0f };
float releaseCoeff{ 0.0f };
float currentGain{ 1.0f };
};
} // namespace FDNReverb

21
Source/DSP/SAPFStage.cpp Normal file
View file

@ -0,0 +1,21 @@
#include "SAPFStage.h"
namespace FDNReverb {
void SAPFStage::prepare(const juce::dsp::ProcessSpec& spec, int delayTargetSamples) {
M = delayTargetSamples;
dl.prepare(spec);
dl.setMaximumDelayInSamples(M + 4);
dl.setDelay(static_cast<float>(M));
}
float SAPFStage::tick(float x) noexcept {
float d = dl.popSample(0);
float w = x + gain * d;
dl.pushSample(0, w);
return d - gain * w;
}
void SAPFStage::reset() noexcept { dl.reset(); }
} // namespace FDNReverb

19
Source/DSP/SAPFStage.h Normal file
View file

@ -0,0 +1,19 @@
#pragma once
#include <JuceHeader.h>
namespace FDNReverb {
class SAPFStage {
public:
void prepare(const juce::dsp::ProcessSpec& spec, int delayTargetSamples);
void setGain(float g) noexcept { gain = juce::jlimit(0.3f, 0.72f, g); }
float tick(float x) noexcept;
void reset() noexcept;
private:
juce::dsp::DelayLine<float, juce::dsp::DelayLineInterpolationTypes::Thiran> dl;
float gain{ 0.618f };
int M{ 0 };
};
} // namespace FDNReverb

197
Source/DSP/Saturator.h Normal file
View file

@ -0,0 +1,197 @@
#pragma once
#include <cmath>
#include <algorithm>
namespace FDNReverb {
enum class SaturationMode {
Warm = 0,
Tape = 1,
Tube = 2,
Hard = 3
};
class Saturator {
public:
Saturator() = default;
void reset() noexcept {
prevInput = 0.0f;
switch (currentMode) {
case SaturationMode::Warm: prevF = 1.0f; break;
case SaturationMode::Tape: prevF = 0.0f; break;
case SaturationMode::Tube: prevF = 1.0f; break;
case SaturationMode::Hard: prevF = 0.0f; break;
}
}
void setMode(SaturationMode mode) noexcept {
if (mode != currentMode) {
currentMode = mode;
reset();
}
}
void setMode(int modeIndex) noexcept {
setMode(static_cast<SaturationMode>(std::clamp(modeIndex, 0, 3)));
}
// -------------------------------------------------------------------------
// * Step B fix: only the drive curve changed; ADAA structure fully preserved
// -------------------------------------------------------------------------
// drive = 1 + amount^3 x 1.0 ( maximum 2.0) -> 1 + amount^2 x 2.5 ( maximum 3.5)
//
// amount | old drive | new drive | effect
// -------|----------|----------|--------------------
// 0.30 | 1.027 | 1.225 | + about 7.5dB stronger
// 0.50 | 1.125 | 1.625 | + about 3.2dB stronger
// 0.70 | 1.343 | 2.225 | + about 4.4dB stronger
// 1.00 | 2.000 | 3.500 | + about 4.9dB stronger
//
// -> plugin 24 harmonics visualization
// -------------------------------------------------------------------------
void setAmount(float amount) noexcept {
amount = std::clamp(amount, 0.0f, 1.0f);
currentAmount = amount;
// * Step B: amount^2 x 2.5 stronger
drive = 1.0f + amount * amount * 2.5f;
wetMix = amount * amount * 0.7f;
dryMix = 1.0f - amount * 0.25f;
}
inline float processSample(float input) noexcept {
if (currentAmount < 1e-4f) return input;
const float dryInput = input;
const float driven = input * drive;
float saturated = 0.0f;
switch (currentMode) {
case SaturationMode::Warm: saturated = processWarm(driven); break;
case SaturationMode::Tape: saturated = processTape(driven); break;
case SaturationMode::Tube: saturated = processTube(driven); break;
case SaturationMode::Hard: saturated = processHard(driven); break;
}
saturated /= drive;
return dryInput * dryMix + saturated * wetMix;
}
private:
// --- Warm: Vicanek x/sqrt(1+x^2) + ADAA 1 ---
inline float processWarm(float x) noexcept {
const float F_x = std::sqrt(1.0f + x * x);
const float dx = x - prevInput;
float y;
constexpr float kTol = 1e-5f;
if (std::abs(dx) < kTol) {
const float xAvg = (x + prevInput) * 0.5f;
y = xAvg / std::sqrt(1.0f + xAvg * xAvg);
}
else {
y = (F_x - prevF) / dx;
}
prevInput = x;
prevF = F_x;
return y;
}
// --- Tape: Pade x(27+x^2)/(27+9x^2) (ADAA intentional ) ---
inline float processTape(float x) noexcept {
if (x > 3.0f) { prevInput = x; return 1.0f; }
if (x < -3.0f) { prevInput = x; return -1.0f; }
const float xsq = x * x;
prevInput = x;
return x * (27.0f + xsq) / (27.0f + 9.0f * xsq);
}
// -------------------------------------------------------------------------
// Tube: asymmetric ADAA + * Step B: kNeg 1.5 -> 2.0
// -------------------------------------------------------------------------
// positive side : f(x) = x/sqrt(1+x^2) F(x) = sqrt(1+x^2)
// negative side : f(x) = x/sqrt(1+(kNeg.x)^2) F(x) = (1/kNeg^2)sqrt(1+(kNeg.x)^2) + fShift
//
// C^1 : x=0 F_pos(0) = F_neg(0) = 1 fShift design
// F_pos(0) = sqrt1 = 1
// F_neg(0) = (1/kNeg^2).sqrt1 + fShift = 1
// -> fShift = 1 - 1/kNeg^2
//
// kNeg=2.0 case : fShift = 1 - 0.25 = 0.75
//
// kNeg stronger effect :
// 2 -> waveform asymmetric
// -> harmonics (2f, 4f) plugin visualization
// -------------------------------------------------------------------------
inline float processTube(float x) noexcept {
// * Step B: kNeg = 1.5f -> 2.0f
constexpr float kNeg = 2.0f;
constexpr float kNeg2 = kNeg * kNeg; // 4.0f
constexpr float invKneg2 = 1.0f / kNeg2; // 0.25f
constexpr float fShift = 1.0f - invKneg2; // 0.75f
float F_x;
if (x >= 0.0f) {
F_x = std::sqrt(1.0f + x * x);
}
else {
const float kx = kNeg * x;
F_x = invKneg2 * std::sqrt(1.0f + kx * kx) + fShift;
}
const float dx = x - prevInput;
const bool signChanged = (x >= 0.0f) != (prevInput >= 0.0f);
float y;
constexpr float kTol = 1e-5f;
if (std::abs(dx) < kTol || signChanged) {
// input -> directly
if (x >= 0.0f) {
y = x / std::sqrt(1.0f + x * x);
}
else {
const float kx = kNeg * x;
y = x / std::sqrt(1.0f + kx * kx);
}
}
else {
y = (F_x - prevF) / dx;
}
prevInput = x;
prevF = F_x;
return y;
}
// --- Hard: clipping + ADAA 1 ---
inline float processHard(float x) noexcept {
float F_x;
if (x > 1.0f) F_x = x - 0.5f;
else if (x < -1.0f) F_x = -x - 0.5f;
else F_x = x * x * 0.5f;
const float dx = x - prevInput;
float y;
constexpr float kTol = 1e-5f;
if (std::abs(dx) < kTol) {
y = std::clamp(x, -1.0f, 1.0f);
}
else {
y = (F_x - prevF) / dx;
}
prevInput = x;
prevF = F_x;
return y;
}
float prevInput{ 0.0f };
float prevF{ 1.0f };
SaturationMode currentMode{ SaturationMode::Warm };
float currentAmount{ 0.0f };
float drive{ 1.0f };
float wetMix{ 0.0f };
float dryMix{ 1.0f };
};
} // namespace FDNReverb

663
Source/DSP/UniversalEngine.cpp Normal file
View file

@ -0,0 +1,663 @@
#include "UniversalEngine.h"
namespace FDNReverb {
namespace {
static bool isMathPrime(int n) noexcept {
if (n < 2) return false;
if (n == 2) return true;
if (n % 2 == 0) return false;
for (int i = 3; i * i <= n; i += 2)
if (n % i == 0) return false;
return true;
}
static int findNearestUniquePrime(int target,
const std::array<int, 16>& usedPrimes,
int usedCount) noexcept {
target = std::max(target, 2);
for (int offset = 0; offset < 100000; ++offset) {
int hi = target + offset;
if (isMathPrime(hi)) {
bool used = false;
for (int k = 0; k < usedCount; ++k)
if (usedPrimes[k] == hi) { used = true; break; }
if (!used) return hi;
}
int lo = target - offset;
if (offset > 0 && lo >= 2 && isMathPrime(lo)) {
bool used = false;
for (int k = 0; k < usedCount; ++k)
if (usedPrimes[k] == lo) { used = true; break; }
if (!used) return lo;
}
}
return target;
}
} // anonymous namespace
UniversalEngine::UniversalEngine() {
fbVec.fill(0.0f);
constexpr float phi = 1.6180339887f;
for (int i = 0; i < FDN_ORDER; ++i) {
lfos[i].state = 12345u + static_cast<uint32_t>(i) * 9876u;
lfos[i].smoothed = 0.0f;
const float angle = static_cast<float>(i) * phi;
const float frac = angle - std::floor(angle);
lfos[i].rateMultiplier = 0.80f + frac * 0.40f;
// * LFO: noise LFO offset
const float cAngle = static_cast<float>(i + 5) * phi;
chorusLFOs[i].phase = cAngle - std::floor(cAngle);
const float cRateAngle = static_cast<float>(i + 11) * phi;
chorusLFOs[i].rateScale = 0.30f + (cRateAngle - std::floor(cRateAngle)) * 0.50f;
}
}
void UniversalEngine::prepare(double sampleRate, int /*maxBlockSize*/) {
fs = sampleRate;
#if AMBIVALENCE_USE_STAGE2_ABSORPTION
MagnitudeResponseFitter::precomputeInteractionMatrix(sampleRate);
#endif
auto getPow2 = [](size_t s) -> size_t {
size_t p = 1;
while (p < s) p *= 2;
return p;
};
size_t totalMemoryNeeded =
getPow2(static_cast<size_t>(fs * 0.5)) // * preDelay (max 500ms)
+ getPow2(static_cast<size_t>(fs * 1.0))
+ getPow2(static_cast<size_t>(fs * 0.05)) * 4
+ getPow2(static_cast<size_t>(fs * 0.5)) * FDN_ORDER
+ getPow2(static_cast<size_t>(fs * 0.05)) * FDN_ORDER * SERIAL_APF_STAGES;
memoryPool.allocate(totalMemoryNeeded);
int mask = 0;
float* ptr = nullptr;
// * PreDelay (max 500ms)
ptr = memoryPool.requestMemory(static_cast<size_t>(fs * 0.5), mask);
preDelayLine.init(ptr, mask);
ptr = memoryPool.requestMemory(static_cast<size_t>(fs * 1.0), mask);
erDelay.init(ptr, mask);
for (int i = 0; i < 4; ++i) {
ptr = memoryPool.requestMemory(static_cast<size_t>(fs * 0.05), mask);
inputDiffusers[i].init(ptr, mask);
}
for (int i = 0; i < FDN_ORDER; ++i) {
ptr = memoryPool.requestMemory(static_cast<size_t>(fs * 0.5), mask);
fdnDelays[i].init(ptr, mask);
for (int s = 0; s < SERIAL_APF_STAGES; ++s) {
ptr = memoryPool.requestMemory(static_cast<size_t>(fs * 0.05), mask);
nestedAllpassDelays[i][s].init(ptr, mask);
}
}
acousticMetrics.prepare(sampleRate, 2000.0f);
currentERTapCount = 0;
currentERDelaySamples.fill(0.0f);
currentERGains.fill(0.0f);
outputLimiter.prepare(sampleRate);
outputEQ.prepare(sampleRate);
duckingAttackCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * 0.010f));
duckingReleaseCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * 0.200f));
duckingEnvelope = 0.0f;
// * DC coefficient : fc ~ 5Hz 1HPF
dcBlockerCoeff = 1.0f - (6.28318530718f * 5.0f / static_cast<float>(fs));
dcX1.fill(0.0f);
dcY1.fill(0.0f);
// * Soft-knee: RMS envelope coefficient (~3ms)
fdnRmsEnv.fill(0.0f);
rmsCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * 0.003f));
reset();
}
void UniversalEngine::reset() {
memoryPool.clear();
fbVec.fill(0.0f);
#if AMBIVALENCE_USE_STAGE2_ABSORPTION
for (auto& lineFilters : absorptionFiltersS2)
for (auto& f : lineFilters) f.reset();
#else
for (auto& f : absorptionFilters) f.reset();
#endif
acousticMetrics.reset();
saturatorL.reset();
saturatorR.reset();
outputLimiter.reset();
outputEQ.reset();
duckingEnvelope = 0.0f;
dcX1.fill(0.0f);
dcY1.fill(0.0f);
fdnRmsEnv.fill(0.0f);
for (auto& dl : fdnDelays) dl.resetState(); // * Thiran allpass state
for (auto& lfo : lfos) lfo.smoothed = 0.0f;
}
void UniversalEngine::setParams(const DSPParams& p) {
activeParams = p;
switch (p.algorithmIndex) {
case 0: case 1: currentTopology = ReverbTopology::Room; break;
case 2: case 3: currentTopology = ReverbTopology::Hall; break;
case 4: currentTopology = ReverbTopology::Plate; break;
case 5: currentTopology = ReverbTopology::Spring; break;
case 6: currentTopology = ReverbTopology::Goldfoil; break;
}
const float attMs = juce::jmax(0.1f, p.duckingAttackMs);
const float relMs = juce::jmax(0.1f, p.duckingRelMs);
duckingAttackCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * attMs * 0.001f));
duckingReleaseCoeff = 1.0f - std::exp(-1.0f / (static_cast<float>(fs) * relMs * 0.001f));
// * PreDelay: ms -> sample count
preDelaySamples = p.preDelayMs * 0.001f * static_cast<float>(fs);
outputEQ.setLoCutHz(p.loCutHz);
outputEQ.setHiCutHz(p.hiCutHz);
updateTopologyAndRouting();
}
void UniversalEngine::calculatePrimePowerDelays() {
const float fsf = static_cast<float>(fs);
const float sizeCoeff = juce::jlimit(0.5f, 2.0f, activeParams.roomSizeScale + 1.0f);
const float minDelayMs = 15.0f + sizeCoeff * 7.5f;
const float maxDelayMs = 50.0f + sizeCoeff * 75.0f;
const int minDelaySamples = std::max(11, static_cast<int>(minDelayMs * 0.001f * fsf));
const int maxDelaySamples = static_cast<int>(maxDelayMs * 0.001f * fsf);
const float logMin = std::log(static_cast<float>(minDelaySamples));
const float logMax = std::log(static_cast<float>(maxDelaySamples));
std::array<int, FDN_ORDER> usedPrimes;
usedPrimes.fill(0);
for (int i = 0; i < FDN_ORDER; ++i) {
const float t = static_cast<float>(i) / static_cast<float>(FDN_ORDER - 1);
const float logTgt = logMin + t * (logMax - logMin);
const int target = static_cast<int>(std::round(std::exp(logTgt)));
const int prime = findNearestUniquePrime(target, usedPrimes, i);
usedPrimes[i] = prime;
fdnBaseDelaySamples[i] = static_cast<float>(prime);
}
}
void UniversalEngine::updateTopologyAndRouting() {
calculatePrimePowerDelays();
auto& preset = *ALL_PRESETS[activeParams.algorithmIndex];
std::array<float, NUM_BANDS> scaledRT60 = preset.acoustics.rt60;
for (auto& v : scaledRT60) v *= activeParams.decayScale;
// -------------------------------------------------------------------------
// * 2) fix : proMode always Tilt / band apply
// -------------------------------------------------------------------------
// old implementation : if (activeParams.proMode) { ... }
// when ProMode is OFF, the Tilt / band coefficients were not applied,
// so the RT60 graph kept the preset's original curve.
//
// new implementation: always apply; the coefficients default to 1.0f,
// so changing them scales the RT60 graph,
// and reset to 1.0f when loadPresetDefaults() is called.
//
// -------------------------------------------------------------------------
scaledRT60[0] *= activeParams.tiltLow;
scaledRT60[1] *= activeParams.tiltLow;
scaledRT60[2] *= activeParams.tiltLow;
scaledRT60[3] *= activeParams.tiltMid;
scaledRT60[4] *= activeParams.tiltMid;
scaledRT60[5] *= activeParams.tiltMid;
scaledRT60[6] *= activeParams.tiltMid;
scaledRT60[7] *= activeParams.tiltHigh;
scaledRT60[8] *= activeParams.tiltHigh;
scaledRT60[9] *= activeParams.tiltHigh;
for (int b = 0; b < NUM_BANDS; ++b)
scaledRT60[b] *= activeParams.rtBands[b];
#if AMBIVALENCE_USE_STAGE2_ABSORPTION
std::array<float, NUM_BANDS> targetDbAccum;
targetDbAccum.fill(0.0f);
for (int i = 0; i < FDN_ORDER; ++i) {
auto s2 = MagnitudeResponseFitter::designStage2(
static_cast<int>(fdnBaseDelaySamples[i]), fs, scaledRT60,
activeParams.hfDamping, activeParams.lfAbsorption);
for (int b = 0; b < NUM_BANDS; ++b) {
currentAbsorptionCoeffsS2[i][b] = s2.geqStages[b];
targetDbAccum[b] += s2.targetDb[b];
}
}
const float representativeDelay = fdnBaseDelaySamples[FDN_ORDER / 2];
for (int b = 0; b < NUM_BANDS; ++b) {
const float avgTargetDb = targetDbAccum[b] / static_cast<float>(FDN_ORDER);
if (avgTargetDb < -0.001f) {
effectiveRT60[b] = -60.0f * representativeDelay
/ (static_cast<float>(fs) * avgTargetDb);
}
else {
effectiveRT60[b] = scaledRT60[b];
}
effectiveRT60[b] = juce::jlimit(0.05f, 30.0f, effectiveRT60[b]);
}
#else
effectiveRT60 = scaledRT60;
for (int i = 0; i < FDN_ORDER; ++i) {
auto absoStages = FilterDesign::designAbsorption(
static_cast<int>(fdnBaseDelaySamples[i]), fs, scaledRT60,
activeParams.hfDamping, activeParams.lfAbsorption);
currentAbsorptionCoeffs[i] = absoStages[0];
}
#endif
// -------------------------------------------------------------------------
// * EDT fix : band average LF/HF correction
// -------------------------------------------------------------------------
// old implementation : effectiveRT60[4] (500Hz) band use
// -> HF Damping high band below EDT
// -> LF Absorption low band below EDT
//
// new implementation : mid-band band (125Hz~4kHz = band 2~7) average value use
// -> band LF/HF correction influence
// -> (31Hz, 63Hz, 8kHz, 16kHz) ( psychoacoustically EDT
// , value unstable )
// -------------------------------------------------------------------------
float rt60Mid = 0.0f;
for (int b = 2; b <= 7; ++b)
rt60Mid += effectiveRT60[b];
rt60Mid = std::max(0.1f, rt60Mid / 6.0f);
// -------------------------------------------------------------------------
// * metallic sound (1): Decay depends on saturation
// -------------------------------------------------------------------------
// each FDN loop pass runs processMicroSaturation(), and reverberation
// nonlinear distortion accumulates in the reverb and shifts the filter
// response, producing metallic ringing.
//
// policy: not applied below a 2.0 s mid-band RT60 average, scaled between 2.0 s and 6.0 s,
// and fully bypassed above 6.0 s.
// -------------------------------------------------------------------------
microSatBlend = juce::jlimit(0.0f, 1.0f, 1.0f - (rt60Mid - 2.0f) / 4.0f);
// -------------------------------------------------------------------------
// * metallic sound (2): Decay depends on modulation
// -------------------------------------------------------------------------
// longer reverb tails require deeper modulation at the filter peaks.
// as used by Lexicon / Strymon.
//
// * modulation depth (scaled down for short reverbs)
// RT60 <= 1.0 s -> 1.0x (min)
// RT60 = 3.0s -> 2.0x
// RT60 >= 5.0 s -> 3.0x (max)
// -------------------------------------------------------------------------
modDepthScale = 1.0f + juce::jlimit(0.0f, 2.0f, (rt60Mid - 1.0f) * 0.5f);
constexpr float baseDB = 16.0f;
float decayCompDB = 7.0f * std::log10(rt60Mid);
static constexpr std::array<float, 7> algorithmOffsetDB = {
+0.8f, +0.9f, +0.5f, +0.5f, +1.5f, +0.6f, +0.6f
};
float algoOffset = algorithmOffsetDB[juce::jlimit(0, 6, activeParams.algorithmIndex)];
switch (currentTopology) {
case ReverbTopology::Room:
bypassER = false; bypassInputDiffusers = false;
apfGain = 0.3f; diffusionSensitivity = 1.0f;
break;
case ReverbTopology::Hall:
bypassER = false; bypassInputDiffusers = false;
apfGain = 0.618f; diffusionSensitivity = 1.0f;
break;
case ReverbTopology::Plate:
bypassER = true; bypassInputDiffusers = false;
apfGain = 0.7f; diffusionSensitivity = 0.7f;
break;
case ReverbTopology::Spring:
bypassER = true; bypassInputDiffusers = false;
apfGain = 0.5f; diffusionSensitivity = 0.5f;
break;
case ReverbTopology::Goldfoil:
bypassER = true; bypassInputDiffusers = false;
apfGain = 0.75f; diffusionSensitivity = 0.8f;
break;
}
const auto& erPattern = PRESET_ER_PATTERNS[
juce::jlimit(0, 6, activeParams.algorithmIndex)];
currentERTapCount = erPattern.numTaps;
float erSizeScale = 0.5f + activeParams.roomSizeScale;
for (int i = 0; i < erPattern.numTaps; ++i) {
currentERDelaySamples[i] = erPattern.taps[i].delayMs * 0.001f
* static_cast<float>(fs) * erSizeScale;
currentERGains[i] = erPattern.taps[i].gain;
}
if (erPattern.numTaps == 0) bypassER = true;
float edtCoeff = 0.7f;
switch (currentTopology) {
case ReverbTopology::Room: edtCoeff = 0.70f; break;
case ReverbTopology::Hall: edtCoeff = 0.95f; break;
case ReverbTopology::Plate: edtCoeff = 0.60f; break;
case ReverbTopology::Spring: edtCoeff = 0.50f; break;
case ReverbTopology::Goldfoil: edtCoeff = 0.85f; break;
}
theoreticalEDT = rt60Mid * edtCoeff;
float satMultiplier = 1.0f;
switch (currentTopology) {
case ReverbTopology::Room: satMultiplier = 0.90f; break;
case ReverbTopology::Hall: satMultiplier = 0.93f; break;
case ReverbTopology::Plate: satMultiplier = 1.00f; break;
case ReverbTopology::Spring: satMultiplier = 1.05f; break;
case ReverbTopology::Goldfoil: satMultiplier = 1.02f; break;
}
float effectiveSatAmount = juce::jlimit(0.0f, 1.0f,
activeParams.saturation * satMultiplier);
saturatorL.setAmount(effectiveSatAmount);
saturatorR.setAmount(effectiveSatAmount);
saturatorL.setMode(activeParams.satTypeIdx);
saturatorR.setMode(activeParams.satTypeIdx);
lateMakeupGainLinear = juce::Decibels::decibelsToGain(baseDB + decayCompDB + algoOffset);
}
inline void UniversalEngine::fastWalshHadamardTransform(
std::array<float, 16>& v) noexcept
{
for (int h = 1; h < 16; h *= 2) {
for (int i = 0; i < 16; i += h * 2) {
for (int j = i; j < i + h; ++j) {
float x = v[j], y = v[j + h];
v[j] = x + y;
v[j + h] = x - y;
}
}
}
for (int i = 0; i < 16; ++i) v[i] *= 0.25f;
}
inline void UniversalEngine::applySignFlipping(
std::array<float, 16>& v) noexcept
{
static constexpr std::array<float, 16> flip = {
1.f, -1.f, 1.f, -1.f, -1.f, 1.f, -1.f, 1.f,
1.f, 1.f, -1.f, -1.f, -1.f, -1.f, 1.f, 1.f
};
for (int i = 0; i < 16; ++i) v[i] *= flip[i];
}
void UniversalEngine::processBlock(const float* inL, const float* inR,
float* outL, float* outR,
int numSamples) noexcept
{
// * CPU: fs float (processBlock throughout use )
const float fsf = static_cast<float>(fs);
// * modulation : squared curve + coefficient suppress
// modAmount^2 low band gradually , 0.001f entire
// : modAmt=0.5 -> 48smp(1ms) / : modAmt=0.5 -> 12smp(0.25ms)
const float modAmtCurved = activeParams.modAmount * activeParams.modAmount;
const float depthSamples = modAmtCurved * 0.001f * fsf * modDepthScale;
const float wetGain = juce::Decibels::decibelsToGain(activeParams.wetDB);
const float stereoWidth = activeParams.stereoWidth;
const float erLevel = activeParams.erLevel;
const float lateLevel = activeParams.lateLevel;
const bool erSolo = activeParams.erSolo;
const float duckThreshLin = juce::Decibels::decibelsToGain(activeParams.duckingThreshDB);
const float duckAmountDB = activeParams.duckingAmount;
const float effectiveDiffusion = activeParams.diffusion * diffusionSensitivity;
const float diffuserGain = 0.25f + effectiveDiffusion * 0.55f;
const float effectiveApfGain = apfGain * (0.60f + effectiveDiffusion * 0.40f);
const float sideBoost = stereoWidth * 1.5f;
const float erLeakage = (1.0f - stereoWidth) * 0.7f;
// * CPU: apfGainStage loop -> before compute
const float apfGainStage = effectiveApfGain * 0.78f;
// * CPU: freqModScale before compute (16ch)
std::array<float, FDN_ORDER> freqModScales;
constexpr float invFdnM1 = 1.0f / static_cast<float>(FDN_ORDER - 1);
for (int i = 0; i < FDN_ORDER; ++i)
freqModScales[i] = 0.5f + (1.0f - static_cast<float>(i) * invFdnM1) * 1.0f;
// * CPU: input diffuser time before compute
std::array<float, 4> diffuserDelaySmp;
for (int i = 0; i < 4; ++i)
diffuserDelaySmp[i] = (3.0f + i * 2.0f) * 0.001f * fsf;
// * CPU: Allpass before compute (16ch x 3)
constexpr float apfBaseMs[SERIAL_APF_STAGES] = { 1.5f, 2.3f, 3.7f };
constexpr float apfSpreadMs[SERIAL_APF_STAGES] = { 0.30f, 0.37f, 0.47f };
constexpr float apfModFrac[SERIAL_APF_STAGES] = { 0.15f, 0.10f, 0.07f };
const float msToSmp = 0.001f * fsf;
std::array<std::array<float, SERIAL_APF_STAGES>, FDN_ORDER> apfBaseDelaySmp;
for (int i = 0; i < FDN_ORDER; ++i)
for (int s = 0; s < SERIAL_APF_STAGES; ++s)
apfBaseDelaySmp[i][s] = (apfBaseMs[s] + i * apfSpreadMs[s]) * msToSmp;
// * CPU: ER tapGain * 0.5f before compute
std::array<float, MAX_ER_TAPS> erTapGainsHalf;
for (int t = 0; t < currentERTapCount; ++t)
erTapGainsHalf[t] = currentERGains[t] * 0.5f;
// * CPU: soft-knee threshold squared before compute (sqrt avoid )
constexpr float compThresh = 0.35f;
constexpr float compThreshSq = compThresh * compThresh;
std::array<float, FDN_ORDER> lfoCoeffs;
{
constexpr float twoPi = 6.28318530718f;
for (int i = 0; i < FDN_ORDER; ++i) {
const float fc = activeParams.modRate * lfos[i].rateMultiplier;
lfoCoeffs[i] = juce::jlimit(0.0001f, 0.9999f,
1.0f - std::exp(-twoPi * fc / fsf));
// * LFO update
chorusLFOs[i].phaseInc = activeParams.modRate * chorusLFOs[i].rateScale / fsf;
}
}
for (int n = 0; n < numSamples; ++n) {
const float leftIn = inL[n];
const float rightIn = inR[n];
const float midIn = (leftIn + rightIn) * 0.5f;
const float sideIn = (leftIn - rightIn) * 0.5f;
float erOutL = 0.0f, erOutR = 0.0f;
// * PreDelay: dry time
// ERFDN input .
// dry attack after ,
// clarity (D50/C50) significantly above .
preDelayLine.write(midIn);
const float delayedMid = (preDelaySamples > 0.5f)
? preDelayLine.read(preDelaySamples)
: midIn;
const float inputPeak = juce::jmax(std::abs(leftIn), std::abs(rightIn));
const float envCoeff = (inputPeak > duckingEnvelope)
? duckingAttackCoeff : duckingReleaseCoeff;
duckingEnvelope += (inputPeak - duckingEnvelope) * envCoeff;
float duckGainLinear = 1.0f;
if (duckAmountDB > 0.001f && duckingEnvelope > duckThreshLin) {
const float envDB = 20.0f * std::log10(juce::jmax(duckingEnvelope, 1e-6f));
const float overDB = envDB - activeParams.duckingThreshDB;
const float gainRedDB = -juce::jmin(overDB, duckAmountDB);
duckGainLinear = juce::Decibels::decibelsToGain(gainRedDB);
}
float fdnInputMid = delayedMid;
if (!bypassInputDiffusers) {
for (int i = 0; i < 4; ++i) {
float d = inputDiffusers[i].read(diffuserDelaySmp[i]);
float w = fdnInputMid + diffuserGain * d;
inputDiffusers[i].write(w);
fdnInputMid = d - diffuserGain * w;
}
}
if (!bypassER) {
erDelay.write(delayedMid);
float erTotalL = 0.0f, erTotalR = 0.0f;
for (int t = 0; t < currentERTapCount; ++t) {
const float tapValue = erDelay.read(currentERDelaySamples[t]);
const float tapGain = erTapGainsHalf[t];
const float tg = tapValue * tapGain;
const float tgLeak = tg * erLeakage;
if (t % 2 == 0) {
erTotalL += tg;
erTotalR += tgLeak;
}
else {
erTotalR += tg;
erTotalL += tgLeak;
}
}
erOutL = erTotalL;
erOutR = erTotalR;
}
// * ER -> Late: feed the ER output into the FDN input
// the early reflections are wall-surface reflections that seed the Late Reverb,
// making the ER-to-Late transition natural and smooth.
if (!bypassER) {
fdnInputMid += (erOutL + erOutR) * 0.5f * 0.15f;
}
std::array<float, 16> currentFb = fbVec;
fastWalshHadamardTransform(currentFb);
applySignFlipping(currentFb);
float fdnOutL = 0.0f, fdnOutR = 0.0f;
std::array<float, 16> nextFb;
for (int i = 0; i < FDN_ORDER; ++i) {
const float lfoVal = lfos[i].tick(lfoCoeffs[i]);
// * modulation: sine-wave LFO + noise LFO
// noise = random (suppresses metallic ringing)
// chorus = smoothly accumulated (rich tail)
const float chorusVal = chorusLFOs[i].tick();
const float combinedLfo = lfoVal + chorusVal * 0.6f;
// * frequency-dependent modulation: high bands modulate less than low bands
const float freqModScale = freqModScales[i];
const float delaySmp = fdnBaseDelaySamples[i]
+ combinedLfo * depthSamples * freqModScale;
float d = fdnDelays[i].read(delaySmp);
#if AMBIVALENCE_USE_STAGE2_ABSORPTION
for (int s = 0; s < ABSO_STAGES_S2; ++s)
d = absorptionFiltersS2[i][s].tick(d, currentAbsorptionCoeffsS2[i][s]);
#else
d = absorptionFilters[i].tick(d, currentAbsorptionCoeffs[i]);
#endif
// * metallic sound (3): DC blocker (1st-order HPF, fc ~ 5 Hz)
// saturation in the FDN loop absorption filters can
// accumulate DC; blocking it prevents low-band asymmetric distortion.
{
const float dcIn = d;
const float dcOut = dcIn - dcX1[i] + dcBlockerCoeff * dcY1[i];
dcX1[i] = dcIn;
dcY1[i] = dcOut;
d = dcOut;
}
// * soft-knee compression (in the FDN feedback loop)
// an RMS envelope over the threshold triggers compression.
// * CPU: sqrt only runs above threshold (compare on squared values)
{
fdnRmsEnv[i] += (d * d - fdnRmsEnv[i]) * rmsCoeff;
if (fdnRmsEnv[i] > compThreshSq) {
const float env = std::sqrt(fdnRmsEnv[i]);
const float over = env - compThresh;
d *= compThresh / (compThresh + over * 0.65f);
}
}
// * metallic sound (1): Decay depends on saturation
// microSatBlend=1.0 -> applied (into the reverb loop)
// microSatBlend=0.0 -> fully bypassed
if (microSatBlend > 0.001f) {
const float sat = processMicroSaturation(d);
d = d + (sat - d) * microSatBlend;
}
// * 3 nested allpass filters (echo density)
// * CPU: apfGainStage precomputed per block
float apfOut = d;
{
for (int s = 0; s < SERIAL_APF_STAGES; ++s) {
const float apfModDepth = depthSamples * apfModFrac[s];
const float apfDelaySmp = apfBaseDelaySmp[i][s]
+ combinedLfo * apfModDepth * freqModScale;
float apfD = nestedAllpassDelays[i][s].read(apfDelaySmp);
float apfW = apfOut + apfGainStage * apfD;
nestedAllpassDelays[i][s].write(apfW);
apfOut = apfD - apfGainStage * apfW;
}
}
nextFb[i] = apfOut;
const float sideForCh = (i % 2 == 0 ? +sideIn : -sideIn) * sideBoost;
const float fdnInputForThisCh = (fdnInputMid + sideForCh) * 0.25f;
fdnDelays[i].write(fdnInputForThisCh + currentFb[i]);
const float crossLeak = 1.0f - stereoWidth;
if (i % 2 == 0) {
fdnOutL += apfOut;
fdnOutR += apfOut * crossLeak;
}
else {
fdnOutR += apfOut;
fdnOutL += apfOut * crossLeak;
}
}
fdnOutL *= 0.125f;
fdnOutR *= 0.125f;
fbVec = nextFb;
const float erMixL = bypassER ? 0.0f : erOutL * erLevel;
const float erMixR = bypassER ? 0.0f : erOutR * erLevel;
const float lateMixL = fdnOutL * lateMakeupGainLinear * lateLevel;
const float lateMixR = fdnOutR * lateMakeupGainLinear * lateLevel;
acousticMetrics.processSample((lateMixL + lateMixR) * 0.5f);
float satL = saturatorL.processSample(lateMixL);
float satR = saturatorR.processSample(lateMixR);
if (erSolo) { satL = 0.0f; satR = 0.0f; }
float wetL = erMixL + satL;
float wetR = erMixR + satR;
outputEQ.process(wetL, wetR);
const float finalWetGain = wetGain * duckGainLinear;
outL[n] = wetL * finalWetGain;
outR[n] = wetR * finalWetGain;
outputLimiter.process(outL[n], outR[n]);
}
}
} // namespace FDNReverb

174
Source/DSP/UniversalEngine.h Normal file
View file

@ -0,0 +1,174 @@
#pragma once
#include "DelayMemory.h"
#include "BiquadFilters.h"
#include "MagnitudeResponseFitter.h"
#include "AcousticMetrics.h"
#include "Saturator.h"
#include "OutputLimiter.h"
#include "OutputEQ.h"
#include "../PluginParameters.h"
#include <array>
#include <cmath>
#define AMBIVALENCE_USE_STAGE2_ABSORPTION 1
namespace FDNReverb {
enum class ReverbTopology { Room, Hall, Plate, Spring, Goldfoil };
// -----------------------------------------------------------------------------
// BandlimitedNoiseLFO: color noise + 1 IIR LPF
// -----------------------------------------------------------------------------
struct BandlimitedNoiseLFO {
uint32_t state{ 12345u };
float smoothed{ 0.0f };
float rateMultiplier{ 1.0f };
inline float nextNoise() noexcept {
state ^= state << 13;
state ^= state >> 17;
state ^= state << 5;
return static_cast<float>(state) * 2.3283064365386963e-10f * 2.0f - 1.0f;
}
inline float tick(float lpfCoeff) noexcept {
smoothed += (nextNoise() - smoothed) * lpfCoeff;
return smoothed;
}
};
// -----------------------------------------------------------------------------
// ChorusLFO: sine-wave phase (modulation)
// -----------------------------------------------------------------------------
struct ChorusLFO {
float phase{ 0.0f };
float phaseInc{ 0.0f };
float rateScale{ 1.0f }; // per-channel rate coefficient (multiplier)
// * CPU: std::sin() replaced by a parabolic approximation (max error ~0.06%, 5-10x faster)
inline float tick() noexcept {
phase += phaseInc;
if (phase >= 1.0f) phase -= 1.0f;
// Parabolic sine: phase [0,1) -> sin(2pi.phase)
const float x = phase < 0.5f ? phase : phase - 1.0f;
const float para = 16.0f * x * (0.5f - std::abs(x));
return para * (0.775f + 0.225f * std::abs(para));
}
};
class UniversalEngine {
public:
UniversalEngine();
void prepare(double sampleRate, int maxBlockSize);
void reset();
void setParams(const DSPParams& p);
void processBlock(const float* inL, const float* inR,
float* outL, float* outR, int numSamples) noexcept;
std::array<float, NUM_BANDS> getEffectiveRT60() const noexcept { return effectiveRT60; }
float getD50() const noexcept { return acousticMetrics.getD50(); }
float getC50() const noexcept { return acousticMetrics.getC50(); }
float getC80() const noexcept { return acousticMetrics.getC80(); }
float getEDT() const noexcept { return theoreticalEDT; }
const AcousticMetrics& getAcousticMetrics() const noexcept { return acousticMetrics; }
int getERTapCount() const noexcept { return currentERTapCount; }
float getERTapDelaySamples(int index) const noexcept {
return (index >= 0 && index < currentERTapCount) ? currentERDelaySamples[index] : 0.0f;
}
float getERTapGain(int index) const noexcept {
return (index >= 0 && index < currentERTapCount) ? currentERGains[index] : 0.0f;
}
double getSampleRate() const noexcept { return fs; }
bool isERBypassed() const noexcept { return bypassER; }
private:
void updateTopologyAndRouting();
void calculatePrimePowerDelays();
inline void fastWalshHadamardTransform(std::array<float, 16>& v) noexcept;
inline void applySignFlipping(std::array<float, 16>& v) noexcept;
// --- FDN loop saturation ---
inline static float processMicroSaturation(float x) noexcept {
constexpr float kInScale = 0.15f;
constexpr float kOutScale = 1.0f / kInScale;
const float xs = x * kInScale;
if (xs > 3.0f) return kOutScale;
if (xs < -3.0f) return -kOutScale;
const float xsq = xs * xs;
return (xs * (27.0f + xsq) / (27.0f + 9.0f * xsq)) * kOutScale;
}
DelayMemoryPool memoryPool;
double fs{ 48000.0 };
DSPParams activeParams;
ReverbTopology currentTopology{ ReverbTopology::Room };
static constexpr int FDN_ORDER = 16;
static constexpr int SERIAL_APF_STAGES = 3; // * Allpass stages
// * PreDelay (max 500 ms)
LinearDelayLine preDelayLine;
float preDelaySamples{ 0.0f };
LinearDelayLine erDelay;
std::array<float, 16> erTaps;
std::array<LinearDelayLine, 4> inputDiffusers;
std::array<ThiranDelayLine, FDN_ORDER> fdnDelays; // * Thiran allpass interpolation
std::array<std::array<LinearDelayLine, SERIAL_APF_STAGES>, FDN_ORDER> nestedAllpassDelays;
int currentERTapCount{ 0 };
std::array<float, MAX_ER_TAPS> currentERDelaySamples;
std::array<float, MAX_ER_TAPS> currentERGains;
OutputLimiter outputLimiter;
OutputEQ outputEQ; // * Phase 5 added
float duckingEnvelope{ 0.0f };
float duckingAttackCoeff{ 0.0f };
float duckingReleaseCoeff{ 0.0f };
#if AMBIVALENCE_USE_STAGE2_ABSORPTION
std::array<std::array<BiquadState, ABSO_STAGES_S2>, FDN_ORDER> absorptionFiltersS2;
std::array<std::array<BiquadCoeffs, ABSO_STAGES_S2>, FDN_ORDER> currentAbsorptionCoeffsS2;
#else
std::array<BiquadState, FDN_ORDER> absorptionFilters;
std::array<BiquadCoeffs, FDN_ORDER> currentAbsorptionCoeffs;
#endif
std::array<BandlimitedNoiseLFO, FDN_ORDER> lfos;
std::array<ChorusLFO, FDN_ORDER> chorusLFOs; // * modulation
std::array<float, FDN_ORDER> fdnBaseDelaySamples;
std::array<float, FDN_ORDER> fbVec;
float apfGain{ 0.618f };
bool bypassER{ false };
bool bypassInputDiffusers{ false }; // * new: default false
float lateMixScale{ 1.0f };
float lateMakeupGainLinear{ 1.0f };
// * Phase 5 addition: Diffusion
float diffusionSensitivity{ 1.0f };
// * metallic sound: DecayTime depends on parameters
float microSatBlend{ 1.0f }; // FDN loop saturation blend (0 = bypass, 1 = full)
float modDepthScale{ 1.0f }; // modulation depth scale (increases with Decay time)
// * DC: prevent DC accumulation in the FDN loop
std::array<float, FDN_ORDER> dcX1;
std::array<float, FDN_ORDER> dcY1;
float dcBlockerCoeff{ 0.999f };
// * soft-knee compression: in the FDN feedback loop
std::array<float, FDN_ORDER> fdnRmsEnv;
float rmsCoeff{ 0.002f };
std::array<float, NUM_BANDS> effectiveRT60;
float theoreticalEDT{ 0.0f };
AcousticMetrics acousticMetrics;
Saturator saturatorL;
Saturator saturatorR;
};
} // namespace FDNReverb