Audacity 3.2.0
EBUR128.cpp
Go to the documentation of this file.
1/**********************************************************************
2
3Audacity: A Digital Audio Editor
4
5EBUR128.cpp
6
7Max Maisel
8
9***********************************************************************/
10
11#include "EBUR128.h"
12#include <cstring>
13
14EBUR128::EBUR128(double rate, size_t channels)
15 : mChannelCount(channels)
16 , mRate(rate)
17{
18 mBlockSize = ceil(0.4 * mRate); // 400 ms blocks
19 mBlockOverlap = ceil(0.1 * mRate); // 100 ms overlap
23 for(size_t channel = 0; channel < mChannelCount; ++channel)
25}
26
28{
29 mSampleCount = 0;
30 mBlockRingPos = 0;
32 memset(mLoudnessHist.get(), 0, HIST_BIN_COUNT*sizeof(long int));
33 for(size_t channel = 0; channel < mChannelCount; ++channel)
34 {
35 mWeightingFilter[channel][0].Reset();
36 mWeightingFilter[channel][1].Reset();
37 }
38}
39
40// fs: sample rate
41// returns array of two Biquads
42//
43// EBU R128 parameter sampling rate adaption after
44// Mansbridge, Stuart, Saoirse Finn, and Joshua D. Reiss.
45// "Implementation and Evaluation of Autonomous Multi-track Fader Control."
46// Paper presented at the 132nd Audio Engineering Society Convention,
47// Budapest, Hungary, 2012."
49{
50 ArrayOf<Biquad> pBiquad(size_t(2), true);
51
52 //
53 // HSF pre filter
54 //
55 double db = 3.999843853973347;
56 double f0 = 1681.974450955533;
57 double Q = 0.7071752369554196;
58 double K = tan(M_PI * f0 / fs);
59
60 double Vh = pow(10.0, db / 20.0);
61 double Vb = pow(Vh, 0.4996667741545416);
62
63 double a0 = 1.0 + K / Q + K * K;
64
65 pBiquad[0].fNumerCoeffs[Biquad::B0] = (Vh + Vb * K / Q + K * K) / a0;
66 pBiquad[0].fNumerCoeffs[Biquad::B1] = 2.0 * (K * K - Vh) / a0;
67 pBiquad[0].fNumerCoeffs[Biquad::B2] = (Vh - Vb * K / Q + K * K) / a0;
68
69 pBiquad[0].fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / a0;
70 pBiquad[0].fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / a0;
71
72 //
73 // HPF weighting filter
74 //
75 f0 = 38.13547087602444;
76 Q = 0.5003270373238773;
77 K = tan(M_PI * f0 / fs);
78
79 pBiquad[1].fNumerCoeffs[Biquad::B0] = 1.0;
80 pBiquad[1].fNumerCoeffs[Biquad::B1] = -2.0;
81 pBiquad[1].fNumerCoeffs[Biquad::B2] = 1.0;
82
83 pBiquad[1].fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
84 pBiquad[1].fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
85
86 return pBiquad;
87}
88
89void EBUR128::ProcessSampleFromChannel(float x_in, size_t channel)
90{
91 double value;
92 value = mWeightingFilter[channel][0].ProcessOne(x_in);
93 value = mWeightingFilter[channel][1].ProcessOne(value);
94 if(channel == 0)
95 mBlockRingBuffer[mBlockRingPos] = value * value;
96 else
97 {
98 // Add the power of additional channels to the power of first channel.
99 // As a result, stereo tracks appear about 3 LUFS louder, as specified.
100 mBlockRingBuffer[mBlockRingPos] += value * value;
101 }
102}
103
105{
108
110 {
111 // A new full block of samples was submitted.
114 }
115 // Close the ring.
117 mBlockRingPos = 0;
118 ++mSampleCount;
119}
120
122{
123 // EBU R128: z_i = mean square without root
124
125 // Calculate Gamma_R from histogram.
126 double sum_v;
127 long int sum_c;
128 HistogramSums(0, sum_v, sum_c);
129
130 // Handle incomplete block if no non-zero block was found.
131 if(sum_c == 0)
132 {
134 HistogramSums(0, sum_v, sum_c);
135 }
136
137 // Histogram values are simplified log(x^2) immediate values
138 // without -0.691 + 10*(...) to safe computing power. This is
139 // possible because they will cancel out anyway.
140 // The -1 in the line below is the -10 LUFS from the EBU R128
141 // specification without the scaling factor of 10.
142 double Gamma_R = log10(sum_v/sum_c) - 1;
143 size_t idx_R = round((Gamma_R - GAMMA_A) * double(HIST_BIN_COUNT) / -GAMMA_A - 1);
144
145 // Apply Gamma_R threshold and calculate gated loudness (extent).
146 HistogramSums(idx_R+1, sum_v, sum_c);
147 if(sum_c == 0)
148 // Silence was processed.
149 return 0;
150 // LUFS is defined as -0.691 dB + 10*log10(sum(channels))
151 return 0.8529037031 * sum_v / sum_c;
152}
153
154void EBUR128::HistogramSums(size_t start_idx, double& sum_v, long int& sum_c)
155{
156 double val;
157 sum_v = 0;
158 sum_c = 0;
159 for(size_t i = start_idx; i < HIST_BIN_COUNT; ++i)
160 {
161 val = -GAMMA_A / double(HIST_BIN_COUNT) * (i+1) + GAMMA_A;
162 sum_v += pow(10, val) * mLoudnessHist[i];
163 sum_c += mLoudnessHist[i];
164 }
165}
166
172void EBUR128::AddBlockToHistogram(size_t validLen)
173{
174 // Reset mBlockRingSize to full state to avoid overflow.
175 // The actual value of mBlockRingSize does not matter
176 // since this is only used to detect if blocks are complete (>= mBlockSize).
178
179 size_t idx;
180 double blockVal = 0;
181 for(size_t i = 0; i < validLen; ++i)
182 blockVal += mBlockRingBuffer[i];
183
184 // Histogram values are simplified log10() immediate values
185 // without -0.691 + 10*(...) to safe computing power. This is
186 // possible because these constant cancel out anyway during the
187 // following processing steps.
188 blockVal = log10(blockVal/double(validLen));
189 // log(blockVal) is within ]-inf, 1]
190 idx = round((blockVal - GAMMA_A) * double(HIST_BIN_COUNT) / -GAMMA_A - 1);
191
192 // idx is within ]-inf, HIST_BIN_COUNT-1], discard indices below 0
193 // as they are below the EBU R128 absolute threshold anyway.
194 if(idx < HIST_BIN_COUNT)
195 ++mLoudnessHist[idx];
196}
#define M_PI
Definition: Distortion.cpp:30
void reinit(Integral count, bool initialize=false)
Definition: MemoryX.h:57
ArrayOf< long int > mLoudnessHist
Definition: EBUR128.h:44
void ProcessSampleFromChannel(float x_in, size_t channel)
Definition: EBUR128.cpp:89
size_t mBlockRingPos
Definition: EBUR128.h:47
void HistogramSums(size_t start_idx, double &sum_v, long int &sum_c)
Definition: EBUR128.cpp:154
void Initialize()
Definition: EBUR128.cpp:27
size_t mChannelCount
Definition: EBUR128.h:51
size_t mBlockOverlap
Definition: EBUR128.h:50
static constexpr double GAMMA_A
EBU R128 absolute threshold.
Definition: EBUR128.h:43
void AddBlockToHistogram(size_t validLen)
Definition: EBUR128.cpp:172
static const size_t HIST_BIN_COUNT
Definition: EBUR128.h:41
size_t mBlockRingSize
Definition: EBUR128.h:48
size_t mBlockSize
Definition: EBUR128.h:49
static ArrayOf< Biquad > CalcWeightingFilter(double fs)
Definition: EBUR128.cpp:48
double mRate
Definition: EBUR128.h:52
double IntegrativeLoudness()
Definition: EBUR128.cpp:121
EBUR128(double rate, size_t channels)
Definition: EBUR128.cpp:14
void NextSample()
Definition: EBUR128.cpp:104
size_t mSampleCount
Definition: EBUR128.h:46
Doubles mBlockRingBuffer
Definition: EBUR128.h:45
ArrayOf< ArrayOf< Biquad > > mWeightingFilter
Definition: EBUR128.h:58
fastfloat_really_inline void round(adjusted_mantissa &am, callback cb) noexcept
Definition: fast_float.h:2512
@ A1
Denominator coefficient indices.
Definition: Biquad.h:29
@ A2
Definition: Biquad.h:29
@ B1
Definition: Biquad.h:27
@ B2
Definition: Biquad.h:27
@ B0
Numerator coefficient indices.
Definition: Biquad.h:27