Implement EBU R128 loudness normalization.

Please note that normalizing to large LUFS values may introduce clipping
of large peaks. However, the recommended target value in EBU R128 of
-23 LUFS should be totally safe but only uses a small amount of the
possible quantization range. For tracks with normal dynamic range,
normalizing to about -10 LUFS should be possible without clipping.

src/effects/Normalize.cpp

@@ -6,6 +6,7 @@
 
   Dominic Mazzoni
   Vaughan Johnson (Preview)
+  Max Maisel (Loudness)
 
 *******************************************************************//**
 
@@ -31,11 +32,12 @@
 
 // Define keys, defaults, minimums, and maximums for the effect parameters
 //
-//     Name       Type     Key                        Def      Min      Max   Scale
+//     Name         Type     Key                        Def      Min      Max   Scale
-Param( Level,     double,  wxT("Level"),               -1.0,    -145.0,  0.0,  1  );
+Param( Level,       double,  wxT("Level"),               -23.0,   -145.0,  0.0,  1  );
-Param( RemoveDC,  bool,    wxT("RemoveDcOffset"),      true,    false,   true, 1  );
+Param( RemoveDC,    bool,    wxT("RemoveDcOffset"),      true,    false,   true, 1  );
-Param( ApplyGain, bool,    wxT("ApplyGain"),           true,    false,   true, 1  );
+Param( ApplyGain,   bool,    wxT("ApplyGain"),           true,    false,   true, 1  );
-Param( StereoInd, bool,    wxT("StereoIndependent"),   false,   false,   true, 1  );
+Param( StereoInd,   bool,    wxT("StereoIndependent"),   false,   false,   true, 1  );
+Param( UseLoudness, bool,    wxT("Use Loudness"),        true,    false,   true, 1  );
 
 BEGIN_EVENT_TABLE(EffectNormalize, wxEvtHandler)
    EVT_CHECKBOX(wxID_ANY, EffectNormalize::OnUpdateUI)
@@ -48,6 +50,7 @@ EffectNormalize::EffectNormalize()
    mDC = DEF_RemoveDC;
    mGain = DEF_ApplyGain;
    mStereoInd = DEF_StereoInd;
+   mUseLoudness = DEF_UseLoudness;
 
    SetLinearEffectFlag(false);
 }
@@ -65,7 +68,7 @@ IdentInterfaceSymbol EffectNormalize::GetSymbol()
 
 wxString EffectNormalize::GetDescription()
 {
-   return _("Sets the peak amplitude of one or more tracks");
+   return _("Sets the peak amplitude or loudness of one or more tracks");
 }
 
 wxString EffectNormalize::ManualPage()
@@ -86,6 +89,7 @@ bool EffectNormalize::DefineParams( ShuttleParams & S ){
    S.SHUTTLE_PARAM( mGain, ApplyGain );
    S.SHUTTLE_PARAM( mDC, RemoveDC );
    S.SHUTTLE_PARAM( mStereoInd, StereoInd );
+   S.SHUTTLE_PARAM( mUseLoudness, UseLoudness );
    return true;
 }
 
@@ -95,6 +99,7 @@ bool EffectNormalize::GetAutomationParameters(CommandParameters & parms)
    parms.Write(KEY_ApplyGain, mGain);
    parms.Write(KEY_RemoveDC, mDC);
    parms.Write(KEY_StereoInd, mStereoInd);
+   parms.Write(KEY_UseLoudness, mUseLoudness);
 
    return true;
 }
@@ -105,11 +110,13 @@ bool EffectNormalize::SetAutomationParameters(CommandParameters & parms)
    ReadAndVerifyBool(ApplyGain);
    ReadAndVerifyBool(RemoveDC);
    ReadAndVerifyBool(StereoInd);
+   ReadAndVerifyBool(UseLoudness);
 
    mLevel = Level;
    mGain = ApplyGain;
    mDC = RemoveDC;
    mStereoInd = StereoInd;
+   mUseLoudness = UseLoudness;
 
    return true;
 }
@@ -145,6 +152,7 @@ bool EffectNormalize::Startup()
          mLevel = -mLevel;
       boolProxy = gPrefs->Read(base + wxT("StereoIndependent"), 0L);
       mStereoInd = (boolProxy == 1);
+      mUseLoudness = false;
 
       SaveUserPreset(GetCurrentSettingsGroup());
 
@@ -163,8 +171,15 @@ bool EffectNormalize::Process()
 
    float ratio;
    if( mGain )
-      // same value used for all tracks
+   {
-      ratio = DB_TO_LINEAR(TrapDouble(mLevel, MIN_Level, MAX_Level));
+      if(mUseLoudness)
+         // LU use 10*log10(...) instead of 20*log10(...)
+         // so multiply level by 2 and use standard DB_TO_LINEAR macro.
+         ratio = DB_TO_LINEAR(TrapDouble(mLevel*2, MIN_Level, MAX_Level));
+      else
+         // same value used for all tracks
+         ratio = DB_TO_LINEAR(TrapDouble(mLevel, MIN_Level, MAX_Level));
+   }
    else
       ratio = 1.0;
 
@@ -234,14 +249,23 @@ bool EffectNormalize::Process()
             // we have a linked stereo track
             // so we need to find it's min, max and offset
             // as they are needed to calc the multiplier for both tracks
+
             track = (WaveTrack *) iter.Next();  // get the next one
             msg =
                topMsg + wxString::Format( _("Analyzing second track of stereo pair: %s"), trackName );
+
             float offset2, extent2;
             bGoodResult = AnalyseTrack(track, msg, progress, offset2, extent2);
             if ( !bGoodResult )
                 break;
-            extent = fmax(extent, extent2);
+
+            if (mUseLoudness)
+               // Loudness: use mean of both tracks.
+               extent = (extent + extent2) / 2;
+            else
+               // Peak: use maximum of both tracks.
+               extent = fmax(extent, extent2);
+
             if( (extent > 0) && mGain )
                mMult = ratio / extent; // we need to use this for both linked tracks
             else
@@ -288,7 +312,7 @@ void EffectNormalize::PopulateOrExchange(ShuttleGui & S)
                                         mDC ? wxT("true") : wxT("false"));
             mDCCheckBox->SetValidator(wxGenericValidator(&mDC));
 
-            S.StartHorizontalLay(wxALIGN_CENTER, false);
+            S.StartHorizontalLay(wxALIGN_LEFT, false);
             {
                mGainCheckBox = S.AddCheckBox(_("Normalize maximum amplitude to"),
                                              mGain ? wxT("true") : wxT("false"));
@@ -306,7 +330,10 @@ void EffectNormalize::PopulateOrExchange(ShuttleGui & S)
             }
             S.EndHorizontalLay();
 
-         
+            mUseLoudnessCheckBox = S.AddCheckBox(_("Use integrative loudness instead of maximum amplitude"),
+                                                 mUseLoudness ? wxT("true") : wxT("false"));
+            mUseLoudnessCheckBox->SetValidator(wxGenericValidator(&mUseLoudness));
+
             mStereoIndCheckBox = S.AddCheckBox(_("Normalize stereo channels independently"),
                                                mStereoInd ? wxT("true") : wxT("false"));
             mStereoIndCheckBox->SetValidator(wxGenericValidator(&mStereoInd));
@@ -350,53 +377,75 @@ bool EffectNormalize::AnalyseTrack(const WaveTrack * track, const wxString &msg,
    bool result = true;
    float min, max;
 
-   if(mGain) {
+   if(mGain)
-      // Since we need complete summary data, we need to block until the OD tasks are done for this track
+   {
-      // TODO: should we restrict the flags to just the relevant block files (for selections)
+      if(mUseLoudness)
-      while (track->GetODFlags()) {
+      {
-         // update the gui
+         CalcEBUR128HPF(track->GetRate());
-         if (ProgressResult::Cancelled == mProgress->Update(
+         CalcEBUR128HSF(track->GetRate());
-            0, _("Waiting for waveform to finish computing...")) )
+         if(mDC)
-            return false;
+         {
-         wxMilliSleep(100);
+            result = AnalyseTrackData(track, msg, progress, ANALYSE_LOUDNESS_DC, offset);
+         }
+         else
+         {
+            result = AnalyseTrackData(track, msg, progress, ANALYSE_LOUDNESS, offset);
+            offset = 0.0;
+         }
+
+         extent = sqrt(mSqSum / mCount.as_double());
       }
+      else
+      {
+         // Since we need complete summary data, we need to block until the OD tasks are done for this track
+         // This is needed for track->GetMinMax
+         // TODO: should we restrict the flags to just the relevant block files (for selections)
+         while (track->GetODFlags()) {
+            // update the gui
+            if (ProgressResult::Cancelled == mProgress->Update(
+               0, _("Waiting for waveform to finish computing...")) )
+               return false;
+            wxMilliSleep(100);
+         }
 
-      // set mMin, mMax.  No progress bar here as it's fast.
+         // set mMin, mMax.  No progress bar here as it's fast.
-      auto pair = track->GetMinMax(mCurT0, mCurT1); // may throw
+         auto pair = track->GetMinMax(mCurT0, mCurT1); // may throw
-      min = pair.first, max = pair.second;
+         min = pair.first, max = pair.second;
 
-   } else {
+         if(mDC)
-      min = -1.0, max = 1.0;   // sensible defaults?
+         {
+            min = -1.0, max = 1.0;   // sensible defaults?
+            result = AnalyseTrackData(track, msg, progress, ANALYSE_DC, offset);
+            min += offset;
+            max += offset;
+         }
+      }
    }
-
+   else if(mDC)
-   if(mDC) {
+   {
-      result = AnalyseDC(track, msg, progress, offset);
+      min = -1.0, max = 1.0;   // sensible defaults?
+      result = AnalyseTrackData(track, msg, progress, ANALYSE_DC, offset);
       min += offset;
       max += offset;
-   } else {
+   }
+   else
+   {
+      min = -1.0, max = 1.0;   // sensible defaults?
       offset = 0.0;
    }
 
-   extent = fmax(fabs(min), fabs(max));
+   if(!mUseLoudness)
+      extent = fmax(fabs(min), fabs(max));
    return result;
 }
 
-//AnalyseDC() takes a track, transforms it to bunch of buffer-blocks,
+//AnalyseTrackData() takes a track, transforms it to bunch of buffer-blocks,
-//and executes AnalyzeData on it...
+//and executes selected AnalyseOperation on it...
-bool EffectNormalize::AnalyseDC(const WaveTrack * track, const wxString &msg,
+bool EffectNormalize::AnalyseTrackData(const WaveTrack * track, const wxString &msg,
-                                double &progress,
+                                double &progress, AnalyseOperation op, float &offset)
-                                float &offset)
 {
    bool rc = true;
 
-   offset = 0.0; // we might just return
-
-   if(!mDC)  // don't do analysis if not doing dc removal
-   {
-      progress += 1.0/double(2*GetNumWaveTracks());
-      return(rc);
-   }
-
    //Transform the marker timepoints to samples
    auto start = track->TimeToLongSamples(mCurT0);
    auto end = track->TimeToLongSamples(mCurT1);
@@ -410,7 +459,8 @@ bool EffectNormalize::AnalyseDC(const WaveTrack * track, const wxString &msg,
    //be shorter than the length of the track being processed.
    Floats buffer{ track->GetMaxBlockSize() };
 
-   mSum = 0.0; // dc offset inits
+   mSum   = 0.0; // dc offset inits
+   mSqSum = 0.0; // rms init
    mCount = 0;
 
    sampleCount blockSamples;
@@ -432,7 +482,12 @@ bool EffectNormalize::AnalyseDC(const WaveTrack * track, const wxString &msg,
       totalSamples += blockSamples;
 
       //Process the buffer.
-      AnalyzeData(buffer.get(), block);
+      if(op == ANALYSE_DC)
+         AnalyseDataDC(buffer.get(), block);
+      else if(op == ANALYSE_LOUDNESS)
+         AnalyseDataLoudness(buffer.get(), block);
+      else if(op == ANALYSE_LOUDNESS_DC)
+         AnalyseDataLoudnessDC(buffer.get(), block);
 
       //Increment s one blockfull of samples
       s += block;
@@ -512,13 +567,46 @@ bool EffectNormalize::ProcessOne(
    return rc;
 }
 
-void EffectNormalize::AnalyzeData(float *buffer, size_t len)
+/// @see AnalyseDataLoudnessDC
+void EffectNormalize::AnalyseDataDC(float *buffer, size_t len)
 {
    for(decltype(len) i = 0; i < len; i++)
       mSum += (double)buffer[i];
    mCount += len;
 }
 
+/// @see AnalyseDataLoudnessDC
+void EffectNormalize::AnalyseDataLoudness(float *buffer, size_t len)
+{
+   float value;
+   for(decltype(len) i = 0; i < len; i++)
+   {
+      value = mR128HSF.ProcessOne(buffer[i]);
+      value = mR128HPF.ProcessOne(value);
+      mSqSum += ((double)value) * ((double)value);
+   }
+   mCount += len;
+}
+
+/// Calculates sample sum (for DC) and EBU R128 weighted square sum
+/// (for loudness). This function has variants which only calculate
+/// sum or square sum for performance improvements if only one of those
+/// values is required.
+/// @see AnalyseDataLoudness
+/// @see AnalyseDataDC
+void EffectNormalize::AnalyseDataLoudnessDC(float *buffer, size_t len)
+{
+   float value;
+   for(decltype(len) i = 0; i < len; i++)
+   {
+      mSum += (double)buffer[i];
+      value = mR128HSF.ProcessOne(buffer[i]);
+      value = mR128HPF.ProcessOne(value);
+      mSqSum += ((double)value) * ((double)value);
+   }
+   mCount += len;
+}
+
 void EffectNormalize::ProcessData(float *buffer, size_t len, float offset)
 {
    for(decltype(len) i = 0; i < len; i++) {
@@ -527,6 +615,46 @@ void EffectNormalize::ProcessData(float *buffer, size_t len, float offset)
    }
 }
 
+// after Juha, https://hydrogenaud.io/index.php/topic,76394.0.html
+void EffectNormalize::CalcEBUR128HPF(float fs)
+{
+   double f0 = 38.13547087602444;
+   double Q  =  0.5003270373238773;
+   double K  = tan(M_PI * f0 / fs);
+
+   mR128HPF.Reset();
+
+   mR128HPF.fNumerCoeffs[Biquad::B0] =  1.0;
+   mR128HPF.fNumerCoeffs[Biquad::B1] = -2.0;
+   mR128HPF.fNumerCoeffs[Biquad::B2] =  1.0;
+
+   mR128HPF.fDenomCoeffs[Biquad::A1] = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
+   mR128HPF.fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
+}
+
+// after Juha, https://hydrogenaud.io/index.php/topic,76394.0.html
+void EffectNormalize::CalcEBUR128HSF(float fs)
+{
+   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;
+
+   mR128HSF.Reset();
+
+   mR128HSF.fNumerCoeffs[Biquad::B0] = (Vh + Vb * K / Q + K * K) / a0;
+   mR128HSF.fNumerCoeffs[Biquad::B1] =       2.0 * (K * K -  Vh) / a0;
+   mR128HSF.fNumerCoeffs[Biquad::B2] = (Vh - Vb * K / Q + K * K) / a0;
+
+   mR128HSF.fDenomCoeffs[Biquad::A1] =   2.0 * (K * K - 1.0) / a0;
+   mR128HSF.fDenomCoeffs[Biquad::A2] = (1.0 - K / Q + K * K) / a0;
+}
+
 void EffectNormalize::OnUpdateUI(wxCommandEvent & WXUNUSED(evt))
 {
    UpdateUI();
@@ -542,10 +670,16 @@ void EffectNormalize::UpdateUI()
    }
    mWarning->SetLabel(wxT(""));
 
+   if (mUseLoudness)
+      mLeveldB->SetLabel(_("LUFS"));
+   else
+      mLeveldB->SetLabel(_("dB"));
+
    // Disallow level stuff if not normalizing
    mLevelTextCtrl->Enable(mGain);
    mLeveldB->Enable(mGain);
    mStereoIndCheckBox->Enable(mGain);
+   mUseLoudnessCheckBox->Enable(mGain);
 
    // Disallow OK/Preview if doing nothing
    EnableApply(mGain || mDC);

src/effects/Normalize.h

@@ -6,6 +6,7 @@
 
   Dominic Mazzoni
   Vaughan Johnson (Preview)
+  Max Maisel (Loudness)
 
 **********************************************************************/
 
@@ -19,6 +20,7 @@
 #include <wx/textctrl.h>
 
 #include "Effect.h"
+#include "Biquad.h"
 
 class ShuttleGui;
 
@@ -58,15 +60,25 @@ public:
 private:
    // EffectNormalize implementation
 
+   enum AnalyseOperation
+   {
+      ANALYSE_DC, ANALYSE_LOUDNESS, ANALYSE_LOUDNESS_DC
+   };
+
    bool ProcessOne(
       WaveTrack * t, const wxString &msg, double& progress, float offset);
    bool AnalyseTrack(const WaveTrack * track, const wxString &msg,
                      double &progress, float &offset, float &extent);
-   void AnalyzeData(float *buffer, size_t len);
+   bool AnalyseTrackData(const WaveTrack * track, const wxString &msg, double &progress,
-   bool AnalyseDC(const WaveTrack * track, const wxString &msg, double &progress,
+                     AnalyseOperation op, float &offset);
-                  float &offset);
+   void AnalyseDataDC(float *buffer, size_t len);
+   void AnalyseDataLoudness(float *buffer, size_t len);
+   void AnalyseDataLoudnessDC(float *buffer, size_t len);
    void ProcessData(float *buffer, size_t len, float offset);
 
+   void CalcEBUR128HPF(float fs);
+   void CalcEBUR128HSF(float fs);
+
    void OnUpdateUI(wxCommandEvent & evt);
    void UpdateUI();
 
@@ -75,11 +87,13 @@ private:
    bool   mGain;
    bool   mDC;
    bool   mStereoInd;
+   bool   mUseLoudness;
 
    double mCurT0;
    double mCurT1;
    float  mMult;
    double mSum;
+   double mSqSum;
    sampleCount    mCount;
 
    wxCheckBox *mGainCheckBox;
@@ -87,9 +101,12 @@ private:
    wxTextCtrl *mLevelTextCtrl;
    wxStaticText *mLeveldB;
    wxStaticText *mWarning;
+   wxCheckBox *mUseLoudnessCheckBox;
    wxCheckBox *mStereoIndCheckBox;
 
    bool mCreating;
+   Biquad mR128HSF;
+   Biquad mR128HPF;
 
    DECLARE_EVENT_TABLE()
 };