Implement EBUR128 integrative loudness class.
This commit is contained in:
parent
a94cda94ae
commit
8555018bd4
|
@ -304,55 +304,6 @@ ArrayOf<Biquad> 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> Biquad::CalcEBUR128WeightingFilter(float fs)
|
||||
{
|
||||
ArrayOf<Biquad> 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)
|
||||
{
|
||||
|
|
|
@ -66,7 +66,6 @@ struct Biquad
|
|||
static ArrayOf<Biquad> CalcButterworthFilter(int order, double fn, double fc, int type);
|
||||
static ArrayOf<Biquad> CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int type);
|
||||
static ArrayOf<Biquad> CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int type);
|
||||
static ArrayOf<Biquad> CalcEBUR128WeightingFilter(float fs);
|
||||
|
||||
static void ComplexDiv (double fNumerR, double fNumerI, double fDenomR, double fDenomI,
|
||||
double* pfQuotientR, double* pfQuotientI);
|
||||
|
|
|
@ -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<Biquad> EBUR128::CalcWeightingFilter(double fs)
|
||||
{
|
||||
ArrayOf<Biquad> 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;
|
||||
}
|
||||
|
||||
|
|
|
@ -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<Biquad> 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<long int> 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<ArrayOf<Biquad>> mWeightingFilter;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
Loading…
Reference in New Issue
Block a user