diff --git a/src/effects/Biquad.cpp b/src/effects/Biquad.cpp index 183a9be69..20c54e8f0 100644 --- a/src/effects/Biquad.cpp +++ b/src/effects/Biquad.cpp @@ -304,55 +304,6 @@ ArrayOf Biquad::CalcChebyshevType2Filter(int order, double fn, double fc return std::move(pBiquad); } -// fs: sample rate -// returns array of two Biquads -// -// EBU R128 parameter sampling rate adaption after -// Mansbridge, Stuart, Saoirse Finn, and Joshua D. Reiss. -// "Implementation and Evaluation of Autonomous Multi-track Fader Control." -// Paper presented at the 132nd Audio Engineering Society Convention, -// Budapest, Hungary, 2012." -ArrayOf Biquad::CalcEBUR128WeightingFilter(float fs) -{ - ArrayOf pBiquad(size_t(2), true); - - // - // HSF pre filter - // - double db = 3.999843853973347; - double f0 = 1681.974450955533; - double Q = 0.7071752369554196; - double K = tan(M_PI * f0 / fs); - - double Vh = pow(10.0, db / 20.0); - double Vb = pow(Vh, 0.4996667741545416); - - double a0 = 1.0 + K / Q + K * K; - - pBiquad[0].fNumerCoeffs[Biquad::B0] = (Vh + Vb * K / Q + K * K) / a0; - pBiquad[0].fNumerCoeffs[Biquad::B1] = 2.0 * (K * K - Vh) / a0; - pBiquad[0].fNumerCoeffs[Biquad::B2] = (Vh - Vb * K / Q + K * K) / a0; - - pBiquad[0].fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / a0; - pBiquad[0].fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / a0; - - // - // HPF weighting filter - // - f0 = 38.13547087602444; - Q = 0.5003270373238773; - K = tan(M_PI * f0 / fs); - - pBiquad[1].fNumerCoeffs[Biquad::B0] = 1.0; - pBiquad[1].fNumerCoeffs[Biquad::B1] = -2.0; - pBiquad[1].fNumerCoeffs[Biquad::B2] = 1.0; - - pBiquad[1].fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K); - pBiquad[1].fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K); - - return std::move(pBiquad); -} - void Biquad::ComplexDiv (double fNumerR, double fNumerI, double fDenomR, double fDenomI, double* pfQuotientR, double* pfQuotientI) { diff --git a/src/effects/Biquad.h b/src/effects/Biquad.h index 29bb8c746..dcb667e2b 100644 --- a/src/effects/Biquad.h +++ b/src/effects/Biquad.h @@ -66,7 +66,6 @@ struct Biquad static ArrayOf CalcButterworthFilter(int order, double fn, double fc, int type); static ArrayOf CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int type); static ArrayOf CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int type); - static ArrayOf CalcEBUR128WeightingFilter(float fs); static void ComplexDiv (double fNumerR, double fNumerI, double fDenomR, double fDenomI, double* pfQuotientR, double* pfQuotientI); diff --git a/src/effects/EBUR128.cpp b/src/effects/EBUR128.cpp index 107c5c5b9..0cf00bfff 100644 --- a/src/effects/EBUR128.cpp +++ b/src/effects/EBUR128.cpp @@ -10,3 +10,171 @@ Max Maisel #include "EBUR128.h" +EBUR128::EBUR128(double rate, size_t channels) + : mChannelCount(channels) + , mRate(rate) +{ + mBlockSize = ceil(0.4 * mRate); // 400 ms blocks + mBlockOverlap = ceil(0.1 * mRate); // 100 ms overlap + mLoudnessHist.reinit(HIST_BIN_COUNT, false); + mBlockRingBuffer.reinit(mBlockSize); + mWeightingFilter.reinit(mChannelCount, false); + for(size_t channel = 0; channel < mChannelCount; ++channel) + mWeightingFilter[channel] = CalcWeightingFilter(mRate); +} + +void EBUR128::Initialize() +{ + mSampleCount = 0; + mBlockRingPos = 0; + mBlockRingSize = 0; + memset(mLoudnessHist.get(), 0, HIST_BIN_COUNT*sizeof(long int)); + for(size_t channel = 0; channel < mChannelCount; ++channel) + { + mWeightingFilter[channel][0].Reset(); + mWeightingFilter[channel][1].Reset(); + } +} + +// fs: sample rate +// returns array of two Biquads +// +// EBU R128 parameter sampling rate adaption after +// Mansbridge, Stuart, Saoirse Finn, and Joshua D. Reiss. +// "Implementation and Evaluation of Autonomous Multi-track Fader Control." +// Paper presented at the 132nd Audio Engineering Society Convention, +// Budapest, Hungary, 2012." +ArrayOf EBUR128::CalcWeightingFilter(double fs) +{ + ArrayOf pBiquad(size_t(2), true); + + // + // HSF pre filter + // + double db = 3.999843853973347; + double f0 = 1681.974450955533; + double Q = 0.7071752369554196; + double K = tan(M_PI * f0 / fs); + + double Vh = pow(10.0, db / 20.0); + double Vb = pow(Vh, 0.4996667741545416); + + double a0 = 1.0 + K / Q + K * K; + + pBiquad[0].fNumerCoeffs[Biquad::B0] = (Vh + Vb * K / Q + K * K) / a0; + pBiquad[0].fNumerCoeffs[Biquad::B1] = 2.0 * (K * K - Vh) / a0; + pBiquad[0].fNumerCoeffs[Biquad::B2] = (Vh - Vb * K / Q + K * K) / a0; + + pBiquad[0].fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / a0; + pBiquad[0].fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / a0; + + // + // HPF weighting filter + // + f0 = 38.13547087602444; + Q = 0.5003270373238773; + K = tan(M_PI * f0 / fs); + + pBiquad[1].fNumerCoeffs[Biquad::B0] = 1.0; + pBiquad[1].fNumerCoeffs[Biquad::B1] = -2.0; + pBiquad[1].fNumerCoeffs[Biquad::B2] = 1.0; + + pBiquad[1].fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K); + pBiquad[1].fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K); + + return std::move(pBiquad); +} + +void EBUR128::ProcessSampleFromChannel(float x_in, size_t channel) +{ + double value; + value = mWeightingFilter[channel][0].ProcessOne(x_in); + value = mWeightingFilter[channel][1].ProcessOne(value); + if(channel == 0) + mBlockRingBuffer[mBlockRingPos] = value * value; + else + { + // Add the power of additional channels to the power of first channel. + // As a result, stereo tracks appear about 3 LUFS louder, as specified. + mBlockRingBuffer[mBlockRingPos] += value * value; + } +} + +void EBUR128::NextSample() +{ + ++mBlockRingPos; + ++mBlockRingSize; + + if(mBlockRingPos % mBlockOverlap == 0) + { + // Process new full block. As incomplete blocks shall be discarded + // according to the EBU R128 specification there is no need for + // some special logic for the last blocks. + if(mBlockRingSize >= mBlockSize) + { + // Reset mBlockRingSize to full state to avoid overflow. + // The actual value of mBlockRingSize does not matter + // since this is only used to detect if blocks are complete (>= mBlockSize). + mBlockRingSize = mBlockSize; + + size_t idx; + double blockVal = 0; + for(size_t i = 0; i < mBlockSize; ++i) + blockVal += mBlockRingBuffer[i]; + + // Histogram values are simplified log10() immediate values + // without -0.691 + 10*(...) to safe computing power. This is + // possible because these constant cancel out anyway during the + // following processing steps. + blockVal = log10(blockVal/double(mBlockSize)); + // log(blockVal) is within ]-inf, 1] + idx = round((blockVal - GAMMA_A) * double(HIST_BIN_COUNT) / -GAMMA_A - 1); + + // idx is within ]-inf, HIST_BIN_COUNT-1], discard indices below 0 + // as they are below the EBU R128 absolute threshold anyway. + if(idx >= 0 && idx < HIST_BIN_COUNT) + ++mLoudnessHist[idx]; + } + } + // Close the ring. + if(mBlockRingPos == mBlockSize) + mBlockRingPos = 0; + ++mSampleCount; +} + +double EBUR128::IntegrativeLoudness() +{ + // EBU R128: z_i = mean square without root + + // Calculate Gamma_R from histogram. + double sum_v = 0; + double val; + long int sum_c = 0; + for(size_t i = 0; i < HIST_BIN_COUNT; ++i) + { + val = -GAMMA_A / double(HIST_BIN_COUNT) * (i+1) + GAMMA_A; + sum_v += pow(10, val) * mLoudnessHist[i]; + sum_c += mLoudnessHist[i]; + } + + // Histogram values are simplified log(x^2) immediate values + // without -0.691 + 10*(...) to safe computing power. This is + // possible because they will cancel out anyway. + // The -1 in the line below is the -10 LUFS from the EBU R128 + // specification without the scaling factor of 10. + double Gamma_R = log10(sum_v/sum_c) - 1; + size_t idx_R = round((Gamma_R - GAMMA_A) * double(HIST_BIN_COUNT) / -GAMMA_A - 1); + + // Apply Gamma_R threshold and calculate gated loudness (extent). + sum_v = 0; + sum_c = 0; + for(size_t i = idx_R+1; i < HIST_BIN_COUNT; ++i) + { + val = -GAMMA_A / double(HIST_BIN_COUNT) * (i+1) + GAMMA_A; + sum_v += pow(10, val) * mLoudnessHist[i]; + sum_c += mLoudnessHist[i]; + } + // LUFS is defined as -0.691 dB + 10*log10(sum(channels)) + return 0.8529037031 * sum_v / sum_c; +} + diff --git a/src/effects/EBUR128.h b/src/effects/EBUR128.h index 037b6537a..883a61d94 100644 --- a/src/effects/EBUR128.h +++ b/src/effects/EBUR128.h @@ -11,14 +11,46 @@ Max Maisel #ifndef __EBUR128_H__ #define __EBUR128_H__ +#include "Biquad.h" #include "MemoryX.h" +#include "SampleFormat.h" /// \brief Implements EBU-R128 loudness measurement. class EBUR128 { public: + EBUR128(double rate, size_t channels); + EBUR128(const EBUR128&) = delete; + EBUR128(EBUR128&&) = delete; + ~EBUR128() = default; + + static ArrayOf CalcWeightingFilter(double fs); + void Initialize(); + void ProcessSampleFromChannel(float x_in, size_t channel); + void NextSample(); + double IntegrativeLoudness(); + inline double IntegrativeLoudnessToLUFS(double loudness) + { return 10 * log10(loudness); } private: + static const size_t HIST_BIN_COUNT = 65536; + /// EBU R128 absolute threshold + static constexpr double GAMMA_A = (-70.0 + 0.691) / 10.0; + ArrayOf mLoudnessHist; + Doubles mBlockRingBuffer; + size_t mSampleCount; + size_t mBlockRingPos; + size_t mBlockRingSize; + size_t mBlockSize; + size_t mBlockOverlap; + size_t mChannelCount; + double mRate; + + /// This is be an array of arrays of the type + /// mWeightingFilter[CHANNEL][FILTER] with + /// CHANNEL = LEFT/RIGHT (0/1) and + /// FILTER = HSF/HPF (0/1) + ArrayOf> mWeightingFilter; }; #endif