Main developements from private repo.

Author: deeuu

src/models/DynamicLoudnessCH2012.cpp

@@ -28,10 +28,11 @@
 #include "../modules/CompressSpectrum.h"
 #include "../modules/WeightSpectrum.h"
 #include "../modules/DoubleRoexBank.h"
+#include "../modules/MultiSourceDoubleRoexBank.h"
+#include "../modules/SpecificPartialLoudnessCHGM2011.h"
 #include "../modules/BinauralInhibitionMG2007.h"
 #include "../modules/InstantaneousLoudness.h"
 #include "../modules/ARAverager.h"
-#include "../modules/PeakFollower.h"
 #include "DynamicLoudnessCH2012.h"
 
 namespace loudness{
@@ -74,11 +75,6 @@ namespace loudness{
         isPresentationDiotic_ = isPresentationDiotic;
     }
 
-    void DynamicLoudnessCH2012::setPeakSTLFollowerUsed(bool isPeakSTLFollowerUsed)
-    {
-        isPeakSTLFollowerUsed_ = isPeakSTLFollowerUsed;
-    }
-
     void DynamicLoudnessCH2012::setBinauralInhibitionUsed(bool isBinauralInhibitionUsed)
     {
         isBinauralInhibitionUsed_ = isBinauralInhibitionUsed;
@@ -124,24 +120,52 @@ namespace loudness{
         pathToFilterCoefs_ = pathToFilterCoefs;
     }
 
+    void DynamicLoudnessCH2012::setPartialLoudnessUsed (bool isPartialLoudnessUsed)
+    {
+        isPartialLoudnessUsed_ = isPartialLoudnessUsed;
+    }
+
+    void DynamicLoudnessCH2012::setWindowSpecGM02 (bool isWindowSpecGM02)
+    {
+        isWindowSpecGM02_ = isWindowSpecGM02;
+    }
+
+    void DynamicLoudnessCH2012::setScalingFactor (Real scalingFactor)
+    {
+        scalingFactor_ = scalingFactor;
+    }
+
+    void DynamicLoudnessCH2012::setAttackTimeSTL (Real attackTimeSTL)
+    {
+        attackTimeSTL_ = attackTimeSTL;
+    }
+
+    void DynamicLoudnessCH2012::setReleaseTimeSTL (Real releaseTimeSTL)
+    {
+        releaseTimeSTL_ = releaseTimeSTL;
+    }
+
     void DynamicLoudnessCH2012::configureModelParameters(const string& setName)
     {
         //common to all
-        setRate(1000);
-        setOuterEarFilter(OME::ANSIS342007_FREEFIELD);
-        setMiddleEarFilter(OME::CHGM2011_MIDDLE_EAR);
-        setSpectrumSampledUniformly(true);
-        setHoppingGoertzelDFTUsed(false);
-        setExcitationPatternInterpolated(false);
-        setInterpolationCubic(true);
-        setSpecificLoudnessOutput(true);
-        setBinauralInhibitionUsed(true);
-        setPresentationDiotic(true);
-        setFirstSampleAtWindowCentre(true);
-        setFilterSpacingInCams(0.1);
-        setCompressionCriterionInCams(0.0);
-        attackTimeSTL_ = 0.016;
-        releaseTimeSTL_ = 0.032;
+        setRate (1000);
+        setOuterEarFilter (OME::ANSIS342007_FREEFIELD);
+        setMiddleEarFilter (OME::CHGM2011_MIDDLE_EAR);
+        setSpectrumSampledUniformly (true);
+        setHoppingGoertzelDFTUsed (false);
+        setExcitationPatternInterpolated (false);
+        setInterpolationCubic (true);
+        setPartialLoudnessUsed (false);
+        setSpecificLoudnessOutput (true);
+        setBinauralInhibitionUsed (true);
+        setPresentationDiotic (true);
+        setFirstSampleAtWindowCentre (true);
+        setFilterSpacingInCams (0.1);
+        setCompressionCriterionInCams (0.0);
+        setWindowSpecGM02 (false);
+        setScalingFactor (1.53e-8);
+        setAttackTimeSTL (0.016);
+        setReleaseTimeSTL (0.032);
         attackTimeLTL_ = 0.1;
         releaseTimeLTL_ = 2.0;
 
@@ -205,7 +229,11 @@ namespace loudness{
         RealVec bandFreqsHz {10, 80, 500, 1250, 2540, 4050, 16001};
 
         //window spec
-        RealVec windowSizeSecs {0.128, 0.064, 0.032, 0.016, 0.008, 0.004};
+        RealVec windowSizeSecs;
+        if (isWindowSpecGM02_)
+            windowSizeSecs = {0.064, 0.032, 0.016, 0.008, 0.004, 0.002};
+        else
+            windowSizeSecs = {0.128, 0.064, 0.032, 0.016, 0.008, 0.004};
         vector<int> windowSizeSamples(6,0);
         //round to nearest sample and force to be even such that centre samples
         //are aligned (using periodic Hann window)
@@ -264,6 +292,8 @@ namespace loudness{
                     (new WeightSpectrum(middleEarFilter_, outerEarFilter_))); 
         }
 
+        int lastSpectrumIdx = modules_.size()-1;
+
         /*
          * Roex filters
          */
@@ -272,18 +302,19 @@ namespace loudness{
         Real doubleRoexBankfactor, instantaneousLoudnessFactor;
         if (isSpecificLoudnessOutput_)
         {
-            doubleRoexBankfactor = 1.53e-8;
+            doubleRoexBankfactor = scalingFactor_;
             instantaneousLoudnessFactor = 1.0;
             LOUDNESS_DEBUG(name_ << ": Excitation pattern will be scaled for specific loudness");
         }
         else
         {
             doubleRoexBankfactor = 1.0;
-            instantaneousLoudnessFactor = 1.53e-8;
+            instantaneousLoudnessFactor = scalingFactor_;
         }
 
-        isBinauralInhibitionUsed_ = isBinauralInhibitionUsed_
-            * (input.getNEars() == 2) * isSpecificLoudnessOutput_;
+        isBinauralInhibitionUsed_ = isBinauralInhibitionUsed_ *
+                                    (input.getNEars() == 2) *
+                                    isSpecificLoudnessOutput_;
         if (isBinauralInhibitionUsed_)
             doubleRoexBankfactor /= 0.75;
 
@@ -294,9 +325,6 @@ namespace loudness{
                                     isExcitationPatternInterpolated_,
                                     isInterpolationCubic_)));
 
-        /*
-         * Binaural inhibition
-         */
         if (isBinauralInhibitionUsed_)
         {
             modules_.push_back(unique_ptr<Module> 
@@ -307,39 +335,68 @@ namespace loudness{
             LOUDNESS_DEBUG(name_ << ": No binaural inhibition.");
         }
         outputModules_["SpecificLoudness"] = modules_.back().get();
-        
-        /*
-        * Instantaneous loudness
-        */   
+
         modules_.push_back(unique_ptr<Module>
                 (new InstantaneousLoudness(instantaneousLoudnessFactor, 
                                            isPresentationDiotic_)));
         outputModules_["InstantaneousLoudness"] = modules_.back().get();
 
-        /*
-         * Short-term loudness
-         */
         modules_.push_back(unique_ptr<Module>
                 (new ARAverager(attackTimeSTL_, releaseTimeSTL_)));
         outputModules_["ShortTermLoudness"] = modules_.back().get();
 
-        /*
-         * Long-term loudness
-         */
         modules_.push_back(unique_ptr<Module>
                 (new ARAverager(attackTimeLTL_, releaseTimeLTL_)));
         outputModules_["LongTermLoudness"] = modules_.back().get();
 
         //configure targets
         configureLinearTargetModuleChain();
 
-        //Option to provide PeakFollower
-        if (isPeakSTLFollowerUsed_)
+        if ((input.getNSources() > 1) && (isPartialLoudnessUsed_))
         {
-            modules_.push_back(unique_ptr<Module> (new PeakFollower(2.0)));
-            outputModules_["PeakShortTermLoudness"] = modules_.back().get();
-            outputModules_["ShortTermLoudness"] -> 
-                addTargetModule (*outputModules_["PeakShortTermLoudness"]);
+            LOUDNESS_DEBUG(name_ 
+                           << ": Setting up modules for partial loudness...");
+
+            modules_.push_back(unique_ptr<Module> 
+                    (new MultiSourceDoubleRoexBank (1.5, 40.2,
+                                    filterSpacingInCams_,
+                                    doubleRoexBankfactor,
+                                    false,
+                                    false)));
+
+            int moduleIdx = modules_.size() - 1;
+
+            // Push spectrum to second excitation transformation stage
+            Module* ptrToWeightedSpectrum = modules_[lastSpectrumIdx].get();
+            ptrToWeightedSpectrum -> addTargetModule (*modules_.back().get());
+
+            outputModules_["MultiSourceExcitation"] = modules_.back().get();
+
+            modules_.push_back(unique_ptr<Module> 
+                    (new SpecificPartialLoudnessCHGM2011()));
+
+            if (isBinauralInhibitionUsed_)
+            {
+                modules_.push_back(unique_ptr<Module>
+                        (new BinauralInhibitionMG2007));
+            }
+
+            outputModules_["SpecificPartialLoudness"] = modules_.back().get();
+
+            modules_.push_back(unique_ptr<Module> 
+                    (new InstantaneousLoudness(1.0, isPresentationDiotic_)));
+            outputModules_["InstantaneousPartialLoudness"] = modules_.back().get();
+
+            modules_.push_back(unique_ptr<Module>
+                (new ARAverager(attackTimeSTL_, releaseTimeSTL_)));
+            outputModules_["ShortTermPartialLoudness"] = modules_.back().get();
+
+            modules_.push_back(unique_ptr<Module>
+                    (new ARAverager(attackTimeLTL_, releaseTimeLTL_)));
+            outputModules_["LongTermPartialLoudness"] = modules_.back().get();
+
+            // configure targets for second (parallel) chain
+            configureLinearTargetModuleChain(moduleIdx);
         }
 
         return 1;

src/models/DynamicLoudnessCH2012.h

@@ -72,15 +72,11 @@ namespace loudness{
      * If you want this setExcitationPatternInterpolated(true); In mode `faster'
      * this is true;
      *
-     * A peak follower can be applied to the short-term loudness using 
-     * setPeakSTLFollowerUsed(true) (default is false).
-     *
      * OUTPUTS:
      *  - "SpecificLoudness"
      *  - "InstantaneousLoudness"
      *  - "ShortTermLoudness"
      *  - "LongTermLoudness"
-     *  - "PeakShortTermLoudness" (optional)
      *
      * REFERENCES:
      *
@@ -131,7 +127,7 @@ namespace loudness{
 
             void setPresentationDiotic(bool isPresentationDiotic);
 
-            void setPeakSTLFollowerUsed(bool isPeakSTLFollowerUsed);
+            void setPartialLoudnessUsed (bool isPartialLoudnessUsed);
 
             void setBinauralInhibitionUsed(bool isBinauralInhibitionUsed);
 
@@ -149,18 +145,27 @@ namespace loudness{
 
             void setPathToFilterCoefs(string pathToFilterCoefs);
 
+            void setWindowSpecGM02 (bool isWindowSpecGM02);
+
+            void setScalingFactor (Real scalingFactor);
+
+            void setAttackTimeSTL (Real attackTimeSTL);
+            void setReleaseTimeSTL (Real releaseTimeSTL);
+
         private:
             virtual bool initializeInternal(const SignalBank &input);
 
             string pathToFilterCoefs_;
             Real filterSpacingInCams_, compressionCriterionInCams_;
             Real attackTimeSTL_, releaseTimeSTL_;
             Real attackTimeLTL_, releaseTimeLTL_;
+            Real scalingFactor_;
             bool isSpectrumSampledUniformly_, isHoppingGoertzelDFTUsed_;
             bool isExcitationPatternInterpolated_;
             bool isInterpolationCubic_, isPresentationDiotic_;
             bool isSpecificLoudnessOutput_, isBinauralInhibitionUsed_;
-            bool isFirstSampleAtWindowCentre_, isPeakSTLFollowerUsed_;
+            bool isFirstSampleAtWindowCentre_;
+            bool isPartialLoudnessUsed_, isWindowSpecGM02_;
             OME::Filter outerEarFilter_, middleEarFilter_;
     }; 
 }

src/models/DynamicLoudnessGM2002.cpp

@@ -36,7 +36,6 @@
 #include "../modules/BinauralInhibitionMG2007.h"
 #include "../modules/InstantaneousLoudness.h"
 #include "../modules/ARAverager.h"
-#include "../modules/PeakFollower.h"
 #include "DynamicLoudnessGM2002.h"
 
 namespace loudness{
@@ -68,11 +67,6 @@ namespace loudness{
     {
         isFirstSampleAtWindowCentre_ = isFirstSampleAtWindowCentre;
     }
-    
-    void DynamicLoudnessGM2002::setPeakSTLFollowerUsed(bool isPeakSTLFollowerUsed)
-    {
-        isPeakSTLFollowerUsed_ = isPeakSTLFollowerUsed;
-    }
 
     void DynamicLoudnessGM2002::setAttackTimeSTL(Real attackTimeSTL)
     {
@@ -191,7 +185,6 @@ namespace loudness{
     {
         //common to all
         setRate(1000);
-        setPeakSTLFollowerUsed(false);
         setOuterEarFilter(OME::ANSIS342007_FREEFIELD);
         setMiddleEarFilter(OME::ANSIS342007_MIDDLE_EAR_HPF);
         setSpectrumSampledUniformly(true);
@@ -208,7 +201,7 @@ namespace loudness{
         setBinauralInhibitionUsed(true);
         setPartialLoudnessUsed(true);
         configureSmoothingTimes("GM2002");
-                
+
         if (setName != "GM2002")
         {
             if (setName == "Faster")
@@ -399,8 +392,10 @@ namespace loudness{
             modules_.push_back(unique_ptr<Module> 
                     (new RoexBankANSIS342007(1.8, 38.9, filterSpacingInCams_)));
         }
+
+
         outputModules_["Excitation"] = modules_.back().get();
-        
+
         /*
          * Specific loudness
          */
@@ -442,15 +437,6 @@ namespace loudness{
         //configure targets
         configureLinearTargetModuleChain();
 
-        //Option to provide PeakFollower
-        if (isPeakSTLFollowerUsed_)
-        {
-            modules_.push_back(unique_ptr<Module> (new PeakFollower(2.0)));
-            outputModules_["PeakShortTermLoudness"] = modules_.back().get();
-            outputModules_["ShortTermLoudness"] -> 
-                addTargetModule (*outputModules_["PeakShortTermLoudness"]);
-        }
-
         // Masking conditions
         if ((input.getNSources() > 1) && (isPartialLoudnessUsed_))
         {

src/models/DynamicLoudnessGM2002.h

@@ -81,9 +81,6 @@ namespace loudness{
      * If you want this setExcitationPatternInterpolated(true). In `Faster'
      * modes, this is true.
      *
-     * A peak follower can be applied to the short-term loudness using 
-     * setPeakSTLFollowerUsed(true) (default is false).
-     *
      * OUTPUTS:
      *  - "Excitation"
      *  - "SpecificLoudness"
@@ -172,8 +169,6 @@ namespace loudness{
 
             void setReleaseTimeLTL(Real releaseTimeLTL);
 
-            void setPeakSTLFollowerUsed(bool isPeakSTLFollowerUsed);
-
             void setOuterEarFilter(const OME::Filter& outerEarFilter);
 
             void setMiddleEarFilter(const OME::Filter& MiddleEarFilter);
@@ -203,7 +198,7 @@ namespace loudness{
             bool isSpectrumSampledUniformly_, isHoppingGoertzelDFTUsed_;
             bool isSpectralResolutionDoubled_, isBinauralInhibitionUsed_;
             bool isSpecificLoudnessANSIS342007_, isFirstSampleAtWindowCentre_;
-            bool isPeakSTLFollowerUsed_, isPartialLoudnessUsed_;
+            bool isPartialLoudnessUsed_;
             string pathToFilterCoefs_;
             OME::Filter outerEarFilter_, middleEarFilter_;
     }; 

src/models/StationaryLoudnessANSIS342007.cpp

@@ -22,7 +22,6 @@
 #include "../modules/MultiSourceRoexBank.h"
 #include "../modules/SpecificPartialLoudnessMGB1997.h"
 #include "../modules/SpecificLoudnessANSIS342007.h"
-#include "../modules/SpecificLoudnessModANSIS342007.h"
 #include "../modules/BinauralInhibitionMG2007.h"
 #include "../modules/InstantaneousLoudness.h"
 #include "StationaryLoudnessANSIS342007.h"
@@ -34,6 +33,7 @@ namespace loudness{
     {
         //Default parameters
         setOuterEarFilter (OME::Filter::ANSIS342007_FREEFIELD);
+        setMiddleEarFilter (OME::Filter::ANSIS342007_MIDDLE_EAR);
         setFilterSpacingInCams (0.1);
         setPresentationDiotic (true);
         setBinauralInhibitionUsed (true);
@@ -64,6 +64,11 @@ namespace loudness{
         outerEarFilter_ = outerEarFilter;
     }
 
+    void StationaryLoudnessANSIS342007::setMiddleEarFilter(const OME::Filter middleEarFilter)
+    {
+        middleEarFilter_ = middleEarFilter;
+    }
+
     void StationaryLoudnessANSIS342007::setFilterSpacingInCams(Real filterSpacingInCams)
     {
         filterSpacingInCams_ = filterSpacingInCams;
@@ -74,18 +79,13 @@ namespace loudness{
         isSpecificLoudnessANSIS342007_ = isSpecificLoudnessANSIS342007;
     }
 
-    void StationaryLoudnessANSIS342007::setAlpha (const RealVec& alpha)
-    {
-        alpha_ = alpha;
-    }
-
     bool StationaryLoudnessANSIS342007::initializeInternal(const SignalBank &input)
     {
         /*
          * Weighting filter
          */
         modules_.push_back(unique_ptr<Module>
-                (new WeightSpectrum(OME::Filter::ANSIS342007_MIDDLE_EAR, outerEarFilter_))); 
+                (new WeightSpectrum(middleEarFilter_, outerEarFilter_))); 
 
         /*
          * Roex filters -> Specific loudness
@@ -96,20 +96,10 @@ namespace loudness{
                 (new RoexBankANSIS342007(1.8, 38.9, filterSpacingInCams_)));
         outputModules_["Excitation"] = modules_.back().get();
 
-        if (!alpha_.empty())
-        {
-            modules_.push_back(unique_ptr<Module>
-                    (new SpecificLoudnessModANSIS342007(
-                        alpha_,
-                        isBinauralInhibitionUsed_)));
-        }
-        else
-        {
-            modules_.push_back(unique_ptr<Module>
-                    (new SpecificLoudnessANSIS342007(
-                        isSpecificLoudnessANSIS342007_,
-                        isBinauralInhibitionUsed_)));
-        }
+        modules_.push_back(unique_ptr<Module>
+                (new SpecificLoudnessANSIS342007(
+                    isSpecificLoudnessANSIS342007_,
+                    isBinauralInhibitionUsed_)));
 
         /*
          * Binaural inhibition
@@ -123,7 +113,7 @@ namespace loudness{
          */   
         modules_.push_back(unique_ptr<Module> 
                 (new InstantaneousLoudness(1.0, isPresentationDiotic_)));
-        outputModules_["InstantaneousLoudness"] = modules_.back().get();
+        outputModules_["Loudness"] = modules_.back().get();
 
         //configure targets
         configureLinearTargetModuleChain();
@@ -136,13 +126,16 @@ namespace loudness{
             // Excitation transformation based on all sources
             modules_.push_back(unique_ptr<Module>
                     (new MultiSourceRoexBank(filterSpacingInCams_)));
-            outputModules_["MultiSourceExcitation"] = modules_.back().get();
+
+            /*
+            modules_.push_back(unique_ptr<Module>
+                (new RoexBankANSIS342007(1.8, 38.9, filterSpacingInCams_)));
+            */
             int moduleIdx = modules_.size() - 1;
 
             // Push spectrum to second excitation transformation stage
             Module* ptrToWeightedSpectrum = modules_[0].get();
-            ptrToWeightedSpectrum -> addTargetModule
-                                     (*outputModules_["MultiSourceExcitation"]);
+            ptrToWeightedSpectrum -> addTargetModule (*modules_.back().get());
 
             // Partial loudness
             modules_.push_back(unique_ptr<Module>
@@ -156,7 +149,7 @@ namespace loudness{
 
             modules_.push_back(unique_ptr<Module> 
                     (new InstantaneousLoudness(1.0, isPresentationDiotic_)));
-            outputModules_["InstantaneousPartialLoudness"] = modules_.back().get();
+            outputModules_["PartialLoudness"] = modules_.back().get();
 
             // configure targets for second (parallel) chain
             configureLinearTargetModuleChain(moduleIdx);

src/models/StationaryLoudnessANSIS342007.h

@@ -100,20 +100,19 @@ namespace loudness{
 
             void setOuterEarFilter(const OME::Filter outerEarFilter);
 
+            void setMiddleEarFilter(const OME::Filter middleEarFilter);
+
             void setFilterSpacingInCams(Real filterSpacingInCams);
 
             void setSpecificLoudnessANSIS342007(bool isSpecificLoudnessANSIS342007_);
 
-            void setAlpha(const RealVec& alpha);
-            
         private:
             virtual bool initializeInternal(const SignalBank &input);
 
             Real filterSpacingInCams_;
             bool isPresentationDiotic_, isPartialLoudnessUsed_, isBinauralInhibitionUsed_;
             bool isSpecificLoudnessANSIS342007_;
-            OME::Filter outerEarFilter_;
-            RealVec alpha_;
+            OME::Filter outerEarFilter_, middleEarFilter_;
     }; 
 }
 

src/models/StationaryLoudnessCHGM2011.cpp

@@ -19,6 +19,8 @@
 
 #include "../modules/WeightSpectrum.h"
 #include "../modules/DoubleRoexBank.h"
+#include "../modules/MultiSourceDoubleRoexBank.h"
+#include "../modules/SpecificPartialLoudnessCHGM2011.h"
 #include "../modules/BinauralInhibitionMG2007.h"
 #include "../modules/InstantaneousLoudness.h"
 #include "StationaryLoudnessCHGM2011.h"
@@ -29,11 +31,12 @@ namespace loudness{
         Model("StationaryLoudnessCHGM2011", false)
     {
         //Default parameters
-        setOuterEarFilter(OME::Filter::ANSIS342007_FREEFIELD);
-        setFilterSpacingInCams(0.1);
-        setPresentationDiotic(true);
-        setBinauralInhibitionUsed(true);
-        setSpecificLoudnessOutput(true);
+        setOuterEarFilter (OME::Filter::ANSIS342007_FREEFIELD);
+        setFilterSpacingInCams (0.1);
+        setPresentationDiotic (true);
+        setBinauralInhibitionUsed (true);
+        setSpecificLoudnessOutput (true);
+        setPartialLoudnessUsed (true);
     }
 
     StationaryLoudnessCHGM2011::~StationaryLoudnessCHGM2011()
@@ -44,6 +47,11 @@ namespace loudness{
         isPresentationDiotic_ = isPresentationDiotic;
     }
 
+    void StationaryLoudnessCHGM2011::setPartialLoudnessUsed(bool isPartialLoudnessUsed)
+    {
+        isPartialLoudnessUsed_ = isPartialLoudnessUsed;
+    }
+
     void StationaryLoudnessCHGM2011::setBinauralInhibitionUsed(bool isBinauralInhibitionUsed)
     {
         isBinauralInhibitionUsed_ = isBinauralInhibitionUsed;
@@ -70,7 +78,8 @@ namespace loudness{
          * Weighting filter
          */
         modules_.push_back(unique_ptr<Module>
-                (new WeightSpectrum(OME::CHGM2011_MIDDLE_EAR, outerEarFilter_))); 
+                (new WeightSpectrum (OME::CHGM2011_MIDDLE_EAR,
+                                     outerEarFilter_))); 
 
         // Set up scaling factors depending on output config
         Real doubleRoexBankfactor, instantaneousLoudnessFactor;
@@ -86,7 +95,10 @@ namespace loudness{
             instantaneousLoudnessFactor = 1.53e-8;
         }
 
-        isBinauralInhibitionUsed_ = isBinauralInhibitionUsed_ * (input.getNEars() == 2) * isSpecificLoudnessOutput_;
+        isBinauralInhibitionUsed_ = isBinauralInhibitionUsed_ *
+                                    (input.getNEars() == 2) *
+                                    isSpecificLoudnessOutput_;
+
         if (isBinauralInhibitionUsed_)
             doubleRoexBankfactor /= 0.75;
 
@@ -116,11 +128,54 @@ namespace loudness{
         modules_.push_back(unique_ptr<Module>
                 (new InstantaneousLoudness(instantaneousLoudnessFactor, 
                                            isPresentationDiotic_)));
-        outputModules_["InstantaneousLoudness"] = modules_.back().get();
+        outputModules_["Loudness"] = modules_.back().get();
 
         //configure targets
         configureLinearTargetModuleChain();
 
+        // Masking conditions
+        if ((input.getNSources() > 1) && (isPartialLoudnessUsed_))
+        {
+            LOUDNESS_DEBUG(name_ 
+                           << ": Setting up modules for partial loudness...");
+
+            modules_.push_back(unique_ptr<Module>
+                    (new MultiSourceDoubleRoexBank (1.5, 40.2,
+                                    filterSpacingInCams_,
+                                    doubleRoexBankfactor,
+                                    false,
+                                    false)));
+
+            int moduleIdx = modules_.size() - 1;
+            // Push spectrum to second excitation transformation stage
+            Module* ptrToWeightedSpectrum = modules_[0].get();
+            ptrToWeightedSpectrum -> addTargetModule (*modules_.back().get());
+
+            modules_.push_back(unique_ptr<Module> 
+                    (new SpecificPartialLoudnessCHGM2011()));
+
+            /*
+            int moduleIdx = modules_.size() - 1;
+            Module* ptrToExcitation = modules_[1].get();
+            ptrToExcitation -> addTargetModule (*modules_.back().get());
+            */
+
+            if (isBinauralInhibitionUsed_)
+            {
+                modules_.push_back(unique_ptr<Module>
+                        (new BinauralInhibitionMG2007));
+            }
+
+            outputModules_["SpecificPartialLoudness"] = modules_.back().get();
+
+            modules_.push_back(unique_ptr<Module> 
+                    (new InstantaneousLoudness(1.0, isPresentationDiotic_)));
+            outputModules_["PartialLoudness"] = modules_.back().get();
+
+            // configure targets for second (parallel) chain
+            configureLinearTargetModuleChain(moduleIdx);
+        }
+
         return 1;
     }
 }

src/models/StationaryLoudnessCHGM2011.h

@@ -73,22 +73,24 @@ namespace loudness{
             StationaryLoudnessCHGM2011();
             virtual ~StationaryLoudnessCHGM2011();
 
-            void setPresentationDiotic(bool isPresentationDiotic);
+            void setPresentationDiotic (bool isPresentationDiotic);
 
-            void setBinauralInhibitionUsed(bool isBinauralInhibitionUsed);
+            void setBinauralInhibitionUsed (bool isBinauralInhibitionUsed);
 
-            void setOuterEarFilter(const OME::Filter outerEarFilter);
+            void setOuterEarFilter (const OME::Filter outerEarFilter);
 
-            void setFilterSpacingInCams(Real filterSpacingInCams);
+            void setFilterSpacingInCams (Real filterSpacingInCams);
 
-            void setSpecificLoudnessOutput(bool isSpecificLoudnessOutput);
+            void setSpecificLoudnessOutput (bool isSpecificLoudnessOutput);
+
+            void setPartialLoudnessUsed (bool setPartialLoudnessUsed);
 
         private:
             virtual bool initializeInternal(const SignalBank &input);
 
             Real filterSpacingInCams_;
             bool isPresentationDiotic_, isBinauralInhibitionUsed_;
-            bool isSpecificLoudnessOutput_;
+            bool isSpecificLoudnessOutput_, isPartialLoudnessUsed_;
             OME::Filter outerEarFilter_;
     }; 
 }

src/models/StationaryLoudnessDIN456311991.cpp

@@ -58,7 +58,7 @@ namespace loudness{
          */   
         modules_.push_back(unique_ptr<Module> 
                 (new InstantaneousLoudnessDIN456311991));
-        outputModules_["InstantaneousLoudness"] = modules_.back().get();
+        outputModules_["Loudness"] = modules_.back().get();
 
         //configure targets
         configureLinearTargetModuleChain();

src/modules/ARAverager.h

@@ -1,6 +1,4 @@
 /*
- * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
- *
  * This file is part of Loudness
  *
  * Loudness is free software: you can redistribute it and/or modify
@@ -52,7 +50,7 @@ namespace loudness{
 
         /** Constructs an ARAverager with attackTime and releaseTime in
          * seconds */
-        ARAverager(Real attackTime, Real releaseTime);
+        ARAverager(Real attackTime=0.1, Real releaseTime=0.1);
 
         virtual ~ARAverager();
 

src/modules/Biquad.cpp

@@ -34,10 +34,37 @@ namespace loudness{
         setACoefs(aCoefs);
     }
 
+    Biquad::Biquad(const string &type) :
+        Module("Biquad"),
+        Filter(2),
+        type_ (type)
+    {}
+
     Biquad::~Biquad() {}
 
     bool Biquad::initializeInternal(const SignalBank &input)
     {
+        if (type_ == "RLB")
+        {
+            RealVec bCoefs = {1.0, -2.0, 1.0};
+            RealVec aCoefs = {1.0, -1.99004745483398, 0.99007225036621};
+            setBCoefs (bCoefs);
+            setACoefs (aCoefs);
+            setCoefficientFs (48000);
+        }
+        else if (type_ == "prefilter")
+        {
+            RealVec bCoefs = {1.53512485958697,
+                              -2.69169618940638, 
+                              1.19839281085285};
+            RealVec aCoefs = {1.0,
+                              -1.69065929318241, 
+                              0.73248077421585};
+            setBCoefs (bCoefs);
+            setACoefs (aCoefs);
+            setCoefficientFs (48000);
+        }
+
         LOUDNESS_ASSERT( bCoefs_.size() == 3 &&
                 aCoefs_.size() == 3,
                 "Filter coefficients do not satisfy filter order");
@@ -72,10 +99,19 @@ namespace loudness{
             bCoefs_[2] = (Vl*omegaSqrd - (Vb*omega/Q) + Vh) / denom;
         }
 
+        /*
+        std::cout << "A" << std::endl;
+        for (int i = 0; i < 3; ++i)
+            std::cout << aCoefs_[i] << std::endl;
+        std::cout << "B" << std::endl;
+        for (int i = 0; i < 3; ++i)
+            std::cout << bCoefs_[i] << std::endl;
+        */
+
         delayLine_.initialize(input.getNSources(),
                               input.getNEars(),
                               input.getNChannels(),
-                              2 * order_,
+                              2,
                               input.getFs());
 
         //output SignalBank
@@ -101,22 +137,16 @@ namespace loudness{
 
                     for (int smp = 0; smp < input.getNSamples(); ++smp)
                     {
-                        //input sample
-                        Real x = inputSignal[smp] * gain_;
-                        
                         //filter
-                        Real y = bCoefs_[0]*x + bCoefs_[1]*z[0] + bCoefs_[2]*z[1]
-                                 - aCoefs_[1]*z[2] - aCoefs_[2]*z[3];
-
-                        //update delay line
-                        z[3] = z[2];
-                        z[2] = y;
-                        z[1] = z[0];
-                        z[0] = x;
-
-                        //output sample
+                        Real x = inputSignal[smp];
+                        Real y = bCoefs_[0] * x + z[0];
+                        z[0] = bCoefs_[1] * x - aCoefs_[1] * y + z[1];
+                        z[1] = bCoefs_[2] * x - aCoefs_[2] * y;
                         outputSignal[smp] = y;
                     }
+
+                    killDenormal (z[0]);
+                    killDenormal (z[1]);
                 }
             }
         }

src/modules/Biquad.h

@@ -36,6 +36,7 @@ namespace loudness{
 
             /** Constructs a Biquad filter with no filter coefficients */
             Biquad();
+            Biquad(const string &type_);
             /** Constructs a Biquad filter using the feedforward and feedback
              * coefficients bCoefs and aCoefs.
              * Both coefficient vectors should each be of size three.
@@ -60,6 +61,7 @@ namespace loudness{
             virtual void processInternal(){};
             virtual void resetInternal();
 
+            std::string type_;
             Real coefficientFs_ = 0;
     };
 }

src/modules/Butter.cpp

@@ -108,7 +108,7 @@ namespace loudness{
                             for (int smp = 0; smp < input.getNSamples(); smp++)
                             {
                                 //input sample
-                                Real x = inputSignal[smp] * gain_;
+                                Real x = inputSignal[smp];
 
                                 //filter
                                 Real y = bCoefs_[0]*(x-z[2]) + bCoefs_[2]*
@@ -126,6 +126,9 @@ namespace loudness{
                                 //output sample
                                 outputSignal[smp] = y;
                             }
+
+                            for (int i = 3; i < 6; ++i)
+                                killDenormal (z[i]);
                         }
                     }
                 }

src/modules/EMA.cpp

@@ -0,0 +1,80 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#include "EMA.h"
+#include "../support/AuditoryTools.h"
+
+namespace loudness{
+
+    EMA::EMA (Real timeConstant, bool asEffectiveAverageTime) :
+        Module("EMA"),
+        timeConstant_ (timeConstant),
+        asEffectiveAverageTime_ (asEffectiveAverageTime)
+    {
+        LOUDNESS_DEBUG(name_ << ": Constructed.");
+    };
+
+    EMA::~EMA()
+    {};
+
+    bool EMA::initializeInternal(const SignalBank &input)
+    {
+        //filter coefficient
+        if (asEffectiveAverageTime_)
+            coef_ = 2.0 / (round (timeConstant_ * input.getFrameRate()) + 1);
+        else
+            coef_ = 1 - exp(-1.0 / (input.getFrameRate() * timeConstant_));
+
+        //output SignalBank
+        output_.initialize(input);
+
+        return 1;
+    }
+
+    void EMA::processInternal(const SignalBank &input)
+    {       
+        int lastSampleIdx = input.getNSamples() - 1;
+
+        for (int src = 0; src < input.getNSources(); src++)
+        {
+            for (int ear = 0; ear < input.getNEars(); ++ear)
+            {
+                for (int chn = 0; chn < input.getNChannels(); ++chn)
+                {
+                    const Real* x = input.getSignalReadPointer(src,
+                            ear, chn, 0);
+                    Real* y = output_.getSignalWritePointer(src,
+                            ear, chn, 0);
+                    killDenormal (y[lastSampleIdx]);
+                    Real yPrev = y[lastSampleIdx];
+
+                    for (int smp = 0; smp < input.getNSamples(); ++smp)
+                    {
+                        y[smp] = coef_ * (x[smp] - yPrev) + yPrev;
+                        yPrev = y[smp];
+                    }
+                }
+            }
+        }
+    }
+
+    //output SignalBanks are cleared so not to worry about filter state
+    void EMA::resetInternal()
+    {}
+}

src/modules/EMA.h

@@ -0,0 +1,49 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#ifndef EMA_H
+#define EMA_H
+
+#include "../support/Module.h"
+
+namespace loudness{
+
+    class EMA : public Module
+    {
+    public:
+ 
+        /** Constructs an EMA with attackTime and releaseTime in
+         * seconds */
+        EMA(Real timeConstant, bool asEffectiveAverageTime=false);
+
+        virtual ~EMA();
+
+    private:
+        virtual bool initializeInternal(const SignalBank &input);
+        virtual bool initializeInternal(){return 0;};
+        virtual void processInternal(const SignalBank &input);
+        virtual void processInternal(){};
+        virtual void resetInternal();
+
+        Real timeConstant_, coef_;
+        bool asEffectiveAverageTime_;
+    };
+}
+
+#endif

src/modules/FIR.cpp

@@ -70,8 +70,7 @@ namespace loudness{
 
                     for (int smp = 0; smp < input.getNSamples(); ++smp)
                     {
-                        //input sample
-                        Real x = inputSignal[smp] * gain_;
+                        Real x = inputSignal[smp];
 
                         //output sample
                         outputSignal[smp] = bCoefs_[0] * x + z[0];

src/modules/FixedRoexBank.cpp

@@ -65,10 +65,12 @@ namespace loudness{
         {
             cam = camLo_ + camStep_ * chn;
             fc = camToHertz (cam);
+            output_.setCentreFreq (chn, fc);
+
             erb = centreFreqToCambridgeERB (fc);
             pu = 4.0 * fc / erb;
+
             t2 = 0.35 * (pu / p51_1k);
-            output_.setCentreFreq (chn, fc);
 
             int j = 0;
             while (j < input.getNChannels())

src/modules/ForwardMaskingPO1998.cpp

@@ -0,0 +1,118 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#include "ForwardMaskingPO1998.h"
+
+namespace loudness{
+
+    ForwardMaskingPO1998::ForwardMaskingPO1998 (Real timeConstant1,
+                                                Real timeConstant2,
+                                                Real weight) :
+        Module ("ForwardMaskingPO1998"),
+        Filter (2),
+        timeConstant1_ (timeConstant1),
+        timeConstant2_ (timeConstant2),
+        weight_ (weight)
+    {
+        LOUDNESS_DEBUG(name_ << ": Constructed.");
+    }
+
+    ForwardMaskingPO1998::~ForwardMaskingPO1998()
+    {}
+
+    bool ForwardMaskingPO1998::initializeInternal(const SignalBank &input)
+    {
+        Real weight2 = 1 - weight_;
+        Real alpha1 = std::exp (-1.0 / (timeConstant1_ * input.getFrameRate()));
+        Real alpha2 = std::exp (-1.0 / (timeConstant2_ * input.getFrameRate()));
+
+        // coefs for  n-1
+        RealVec bCoefs = {-(weight2 * alpha2 + weight_ * alpha1)};
+        // coefs for n, n-1 and n-2
+        RealVec aCoefs = {1.0, -alpha1 - alpha2, alpha1 * alpha2};
+
+        aCoefs[0] = 1.0 / ((1 - weight_) * timeConstant1_ + 
+                    weight_ * timeConstant2_);
+        aCoefs[0] /= input.getFrameRate();
+
+        setBCoefs (bCoefs);
+        setACoefs (aCoefs);
+
+        delayLine_.initialize (
+                input.getNSources(),
+                input.getNEars(),
+                input.getNChannels(),
+                2,
+                input.getFs());
+                    
+
+        //output SignalBank
+        output_.initialize(input);
+
+        return 1;
+    }
+
+    void ForwardMaskingPO1998::processInternal(const SignalBank &input)
+    {       
+        for (int src = 0; src < input.getNSources(); ++src)
+        {
+            for(int ear = 0; ear < input.getNEars(); ++ear)
+            {
+                for (int chn = 0; chn < input.getNChannels(); ++chn)
+                {
+                    const Real* inputSignal = input.getSignalReadPointer
+                                              (src, ear, chn);
+                    Real* outputSignal = output_.getSignalWritePointer
+                                         (src, ear, chn);
+                    Real* z = delayLine_.getSignalWritePointer
+                                         (src, ear, chn);
+
+                    for (int smp = 0; smp < input.getNSamples(); ++smp)
+                    {
+                        //input sample
+                        Real x = inputSignal[smp];
+
+                        //output sample
+                        outputSignal[smp] = x + z[0];
+
+                        // b[n-1] ...
+                        z[0] = bCoefs_[0] * x + z[1];
+                        z[1] = 0.0;
+
+                        // a[n-1] ...
+                        z[0] -= aCoefs_[1] * outputSignal[smp];
+
+                        // a[n-2] ...
+                        z[1] -= aCoefs_[2] * outputSignal[smp];
+                        
+                        outputSignal[smp] *= aCoefs_[0];
+                            
+                    }
+
+                    killDenormal (z[0]);
+                    killDenormal (z[1]);
+                }
+            }
+        }
+    }
+
+    //output SignalBanks are cleared so not to worry about filter state
+    void ForwardMaskingPO1998::resetInternal()
+    {}
+}

src/modules/ForwardMaskingPO1998.h

@@ -0,0 +1,48 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#ifndef FORWARDMASKINGPO1998_H
+#define FORWARDMASKINGPO1998_H
+
+#include "../support/Filter.h"
+
+namespace loudness{
+
+    class ForwardMaskingPO1998 : public Module, public Filter
+    {
+        public:
+            ForwardMaskingPO1998(
+                    Real timeConstant1=4e-3,
+                    Real timeConstant2=29e-3,
+                    Real weight=0.0251);
+
+            virtual ~ForwardMaskingPO1998();
+
+        private:
+            virtual bool initializeInternal(const SignalBank &input);
+            virtual bool initializeInternal(){return 0;};
+            virtual void processInternal(const SignalBank &input);
+            virtual void processInternal(){};
+            virtual void resetInternal();
+
+            Real timeConstant1_, timeConstant2_, weight_;
+    };
+}
+
+#endif

src/modules/HoppingGoertzelDFT.cpp

@@ -245,6 +245,12 @@ namespace loudness{
                                 vPrev[k] = v;
                             }
                         }
+
+                        for (int k = binIdxForGoertzels_[w][0]; k < binIdxForGoertzels_[w][1]; ++k)
+                        {
+                            killDenormal (vPrev[k]);
+                            killDenormal (vPrev2[k]);
+                        }
                     }
                 }
             }
@@ -299,7 +305,6 @@ namespace loudness{
             nSamplesUntilTrigger_ = largestWindowSize_ / 2;
         else
             nSamplesUntilTrigger_ = largestWindowSize_;
-        nSamplesUntilTrigger_ = hopSize_;
 
         for (int w = 0; w < nWindows_; ++w)
         {

src/modules/IIR.cpp

@@ -94,7 +94,7 @@ namespace loudness{
                     for (int smp = 0; smp < input.getNSamples(); ++smp)
                     {
                         //input sample
-                        Real x = inputSignal[smp] * gain_;
+                        Real x = inputSignal[smp];
 
                         //output sample
                         outputSignal[smp] = bCoefs_[0] * x + z[0];
@@ -108,6 +108,9 @@ namespace loudness{
                             z[j-1] -= aCoefs_[j] * outputSignal[smp];
                         z[orderMinus1_] -= aCoefs_[order_] * outputSignal[smp];
                     }
+
+                    for (int j = 0; j < order_; ++j)
+                        killDenormal (z[j]);
                 }
             }
         }

src/modules/MultiSourceDoubleRoexBank.cpp

@@ -0,0 +1,266 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#include "MultiSourceDoubleRoexBank.h"
+#include "../support/AuditoryTools.h"
+
+namespace loudness{
+
+    MultiSourceDoubleRoexBank::MultiSourceDoubleRoexBank(Real camLo, 
+            Real camHi,
+            Real camStep,
+            Real scalingFactor,
+            bool isExcitationPatternInterpolated,
+            bool isInterpolationCubic) : 
+        Module("MultiSourceDoubleRoexBank"),
+        camLo_(camLo), 
+        camHi_(camHi),
+        camStep_(camStep),
+        scalingFactor_(scalingFactor),
+        isExcitationPatternInterpolated_(isExcitationPatternInterpolated),
+        isInterpolationCubic_(isInterpolationCubic)
+    {}
+
+    MultiSourceDoubleRoexBank::~MultiSourceDoubleRoexBank() {}
+
+    bool MultiSourceDoubleRoexBank::initializeInternal(const SignalBank &input)
+    {
+
+        if (camStep_ <= 0.1)
+            isExcitationPatternInterpolated_ = false;
+        if (isExcitationPatternInterpolated_)
+        {
+            camLo_ = 1.5;
+            camHi_ = 40.2;
+        }
+
+        //number of roex filters to use
+        nFilters_ = std::floor((camHi_ - camLo_) / camStep_) + 1;
+
+        LOUDNESS_DEBUG(name_
+                << ": interpolation applied: " << isExcitationPatternInterpolated_
+                << " filter spacing in Cams: " << camStep_
+                << " Total number of filters: " << nFilters_);
+
+        //initialise output SignalBank
+        if (isExcitationPatternInterpolated_)
+        {
+            //centre freqs in cams
+            cams_.assign (nFilters_, 0.0);
+
+            //required for log interpolation
+            logExcitation_.assign (nFilters_, 0.0);
+
+            //388 filters to cover [1.5, 40.2] see p. 3
+            output_.initialize (input.getNSources(),
+                                input.getNEars(),
+                                388,
+                                1,
+                                input.getFs());
+            output_.setChannelSpacingInCams (0.1);
+            for (int i = 0; i < 388; ++i)
+                output_.setCentreFreq (i, camToHertz (camLo_ + (i * 0.1)));
+        }
+        else
+        {
+            output_.initialize (
+                    input.getNSources(),
+                    input.getNEars(),
+                    nFilters_,
+                    1,
+                    input.getFs());
+            output_.setChannelSpacingInCams (camStep_);
+        }
+
+        output_.setFrameRate (input.getFrameRate());
+
+        //filter variables
+        wPassive_.resize (nFilters_);
+        wActive_.resize (nFilters_);
+        maxGdB_.resize (nFilters_);
+        thirdGainTerm_.resize (nFilters_);
+
+        //fill the above arrays
+        for (int i = 0; i < nFilters_; ++i)
+        {
+            //filter frequency in Cams
+            Real cam = camLo_ + (i * camStep_);
+            //filter frequency in Hz
+            Real fc = camToHertz (cam);
+
+            if (isExcitationPatternInterpolated_)
+                cams_[i] = cam;
+            else
+                output_.setCentreFreq (i, fc);
+
+            //slopes
+            Real tl = fc / (0.108 * fc + 2.33);
+            Real tu = 15.6;
+            Real pl = fc / (0.027 * fc + 5.44);
+            Real pu = 27.9;
+
+            //precalculate some gain terms
+            maxGdB_[i] = fc / (0.0191 * fc + 1.1);
+            //Third term in Eq. 6
+            thirdGainTerm_[i] = maxGdB_[i] / (1 + exp (0.05 * (100 - maxGdB_[i])));
+
+            //compute the fixed filters
+            int j = 0;
+            while (j < input.getNChannels())
+            {
+                Real pgPassive = 0.0, pgActive = 0.0;
+
+                //normalised deviation
+                Real g = (input.getCentreFreq(j) - fc) / fc;
+
+                //Is g limited to 2 sufficient for the passive filter?
+                if (g <= 2)
+                {
+                    if (g < 0) //lower value 
+                    {
+                        pgPassive = -tl * g; 
+                        pgActive = -pl * g; 
+                    }
+                    else //upper value
+                    {
+                        pgPassive = tu * g;
+                        pgActive = pu * g; 
+                    }
+
+                    //Eq. 4 and Eq. 7
+                    wPassive_[i].push_back ((1 + pgPassive) * exp (-pgPassive)); 
+                    wActive_[i].push_back ((1 + pgActive) * exp (-pgActive)); 
+                }
+                else
+                    break;
+                j++;
+            }
+        }
+        LOUDNESS_DEBUG(name_ << ": Passive and active filters configured.");
+        LOUDNESS_DEBUG(name_ << ": Excitation pattern will be scaled by: " 
+                << scalingFactor_);
+        
+        return 1;
+    }
+
+
+    void MultiSourceDoubleRoexBank::processInternal(const SignalBank &input)
+    {
+        for (int ear = 0; ear < input.getNEars(); ++ear)
+        {
+            // First to passive filtering
+            RealVec gainOut (nFilters_, 0.0);
+            for (int src = 0; src < input.getNSources(); ++src)
+            {
+                const Real* inputSpectrum = input
+                                            .getSingleSampleReadPointer
+                                            (src, ear, 0);
+
+                Real* outputExcitation = output_
+                                         .getSingleSampleWritePointer
+                                         (src, ear, 0);
+
+                for (int i = 0; i < nFilters_; ++i)
+                {
+                    Real excitationLinP = 0.0;
+
+                    //passive filter output
+                    for (uint j = 0; j < wPassive_[i].size(); ++j)
+                        excitationLinP += wPassive_[i][j] * inputSpectrum[j];
+
+                    outputExcitation[i] = excitationLinP;
+
+                    // Accumulate power in each channel from all sources
+                    gainOut[i] += excitationLinP;
+                }
+            }
+
+            // Now gain calculation using total power from all inputs
+            for (int i = 0; i < nFilters_; ++i)
+            {
+                //convert to dB
+                Real excitationLog = powerToDecibels (gainOut[i]);
+
+                //compute gain (Complete Eq. 6 for <= 30)
+                Real gain = maxGdB_[i] - (maxGdB_[i] / 
+                        (1 + exp (-0.05*(excitationLog - (100 - maxGdB_[i]))))) +
+                        thirdGainTerm_[i];
+
+                //check for higher levels
+                if (excitationLog > 30)
+                {
+                    //complete Eq. 6 for > 30
+                    Real excitationLogMinus30 = excitationLog - 30;
+                    gain = gain - 0.003 * excitationLogMinus30 * excitationLogMinus30;
+                }
+
+                //convert to linear gain
+                gainOut[i] = decibelsToPower(gain);
+            }
+            
+            // Now the active filter for each source
+            for (int src = 0; src < input.getNSources(); ++src)
+            {
+                const Real* inputSpectrum = input
+                                            .getSingleSampleReadPointer
+                                            (src, ear, 0);
+
+                Real* outputExcitation = output_.
+                                         getSingleSampleWritePointer
+                                         (src, ear, 0);
+
+                for (int i = 0; i < nFilters_; ++i)
+                {
+                    // from above
+                    Real excitationLinP = outputExcitation[i];
+                    Real excitationLinA = 0.0;
+
+                    //active filter output
+                    for (uint j = 0; j < wActive_[i].size(); ++j)
+                        excitationLinA += wActive_[i][j] * inputSpectrum[j];
+
+                    // Use the gain derived from all inputs
+                    excitationLinA *= gainOut[i];
+
+                    //excitation pattern
+                    Real excitation = scalingFactor_ * 
+                                      (excitationLinP + excitationLinA);
+
+                    if (isExcitationPatternInterpolated_)
+                        logExcitation_[i] = log(excitation + 1e-10);
+                    else
+                        outputExcitation[i] = excitation;
+                }
+
+                //Interpolate to estimate 0.1~Cam res excitation pattern
+                if (isExcitationPatternInterpolated_)
+                {
+                    spline_.set_points (cams_,
+                                        logExcitation_,
+                                        isInterpolationCubic_);
+
+                    for (int i = 0; i < 388; ++i)
+                        outputExcitation[i] = exp (spline_ (camLo_ + i * 0.1));
+                }
+            }
+        }
+    }
+
+    void MultiSourceDoubleRoexBank::resetInternal(){};
+}

src/modules/MultiSourceDoubleRoexBank.h

@@ -0,0 +1,98 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#ifndef MultiSourceDoubleRoexBank_H
+#define MultiSourceDoubleRoexBank_H
+
+#include "../support/Module.h"
+#include "../thirdParty/spline/Spline.h"
+
+namespace loudness{
+
+    /**
+     * @class MultiSourceDoubleRoexBank
+     *
+     * @brief Applies a bank of double roex filters to an input power spectrum.
+     *
+     * This is an implementation of the double roex filter model proposed by
+     * Chen et al (2011). The model makes use of a passive wideband filter and
+     * an active narrowband filter.  The output of the passive filter controls
+     * the active gain of the narrowband filter.
+     *
+     * The centre frequencies of the filters are calculated according to the
+     * Hertz to Cam equation given in ANSI S3.4:2007. Therefore, they will be
+     * slightly different to those computed in Chen et al. (2011) as they use
+     * the traditional equation with rounded coefficients. For example, 40.2
+     * Cams ~= 17197 Hz (ANSI) and 40.2 Cams ~= 17070 Hz (Chen et al.).
+     *
+     * REFERENCES:
+     *
+     * Chen, Z., Hu, G., Glasberg, B. R., & Moore, B. C. J. (2011). A new method
+     * of calculating auditory excitation patterns and loudness for steady
+     * sounds.  Hearing Research, 282(1-2), 204–15.
+     *
+     */
+ 
+    class MultiSourceDoubleRoexBank : public Module
+    {
+
+    public:
+
+        /** Constructs a MultiSourceDoubleRoexBank with double roex filters equally spaced
+         * on the Cam scale.
+         * 
+         * @param camLo Frequency of the first double roex filter in Cams.
+         * @param camHi Frequency of the last double roex filter in Cams.
+         * @param camStep Interval between adjacent filters on the cam Scale.
+         * @param scalingFactor A scalar applied to the output excitation
+         * pattern (default is 1).
+         * @param isExcitationPatternInterpolated Set true to interpolate the
+         * excitation pattern.  
+         * @param isInterpolationCubic Set true for cubic
+         * interpolation, false for linear. Only applies if
+         * isExcitationPatternInterpolated is true;
+         
+         */
+        MultiSourceDoubleRoexBank(Real camLo = 1.5,
+                Real camHi = 40.2,
+                Real camStep = 0.1,
+                Real scalingFactor = 1.0,
+                bool isExcitationPatternInterpolated = false,
+                bool isInterpolationCubic = true);
+
+        virtual ~MultiSourceDoubleRoexBank();
+
+    private:
+
+        virtual bool initializeInternal(const SignalBank &input);
+        virtual bool initializeInternal(){return 0;};
+        virtual void processInternal(const SignalBank &input);
+        virtual void processInternal(){};
+        virtual void resetInternal();
+
+        Real camLo_, camHi_, camStep_, scalingFactor_;
+        bool isExcitationPatternInterpolated_, isInterpolationCubic_;
+        int nFilters_;
+        RealVec maxGdB_, thirdGainTerm_, cams_, logExcitation_;
+        RealVecVec wPassive_, wActive_;
+        spline spline_;
+    };
+}
+
+#endif

src/modules/PeakFollower.cpp

@@ -52,25 +52,28 @@ namespace loudness{
     {       
         int lastSampleIdx = input.getNSamples() - 1;
 
-        for (int ear = 0; ear < input.getNEars(); ++ear)
+        for (int src = 0; src < input.getNSources(); ++src)
         {
-            for (int chn = 0; chn < input.getNChannels(); ++chn)
+            for (int ear = 0; ear < input.getNEars(); ++ear)
             {
-                const Real* x = input.getSignalReadPointer(ear, chn, 0);
+                for (int chn = 0; chn < input.getNChannels(); ++chn)
+                {
+                    const Real* x = input.getSignalReadPointer(src, ear, chn, 0);
 
-                Real* y = output_.getSignalWritePointer(ear, chn, 0);
+                    Real* y = output_.getSignalWritePointer(src, ear, chn, 0);
 
-                Real yPrev = y[lastSampleIdx];
+                    Real yPrev = y[lastSampleIdx];
 
-                for (int smp = 0; smp < input.getNSamples(); ++smp)
-                {
-                    Real absX = abs (x[smp]);
+                    for (int smp = 0; smp < input.getNSamples(); ++smp)
+                    {
+                        Real absX = abs (x[smp]);
 
-                    if (absX >= yPrev)
-                        y[smp] = absX;
-                    else
-                        y[smp] = absX + coef_ * (yPrev - absX);
-                    yPrev = y[smp];
+                        if (absX >= yPrev)
+                            y[smp] = absX;
+                        else
+                            y[smp] = absX + coef_ * (yPrev - absX);
+                        yPrev = y[smp];
+                    }
                 }
             }
         }

src/modules/SpecificPartialLoudnessCHGM2011.cpp

@@ -0,0 +1,99 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#include "SpecificPartialLoudnessCHGM2011.h"
+#include "../support/AuditoryTools.h"
+
+namespace loudness{
+
+    SpecificPartialLoudnessCHGM2011::SpecificPartialLoudnessCHGM2011() :
+        Module("SpecificPartialLoudnessCHGM2011")
+    {
+        LOUDNESS_DEBUG(name_ << ": Constructed.");
+    }
+
+    SpecificPartialLoudnessCHGM2011::~SpecificPartialLoudnessCHGM2011()
+    {};
+
+    bool SpecificPartialLoudnessCHGM2011::initializeInternal(const SignalBank &input)
+    {
+        for (int chn = 0; chn < input.getNChannels(); ++chn)
+        {
+            Real fc = input.getCentreFreq (chn);
+            Real ef = hertzToNormalisedCam (fc);
+            // Table 1 of Chen et al., 2011
+            Real kDB = 6.51 * ef*ef - 1.93;
+            LOUDNESS_DEBUG(fc << ": " << kDB);
+            k_.push_back (std::pow (10, kDB / 10.0));
+        }
+
+        output_.initialize (input);
+
+        return 1;
+    }
+
+    void SpecificPartialLoudnessCHGM2011::processInternal(const SignalBank &input)
+    {       
+        for (int ear = 0; ear < input.getNEars(); ++ear)
+        {
+            /*
+             * Should probably calculate excitation patterns using total power
+             * from all sources? They don't specify this in the paper...
+             * Just go simple linear approach for now.
+             */
+            RealVec eTot (input.getNChannels(), 0.0);
+            for (int src = 0; src < input.getNSources(); ++src)
+            {
+                const Real* inputExcitation = input
+                                              .getSingleSampleReadPointer
+                                              (src, ear, 0);
+
+                for (int chn = 0; chn < input.getNChannels(); ++chn)
+                    eTot[chn] += inputExcitation[chn];
+            }
+
+            // Now do partial loudness calculation
+            for (int src = 0; src < input.getNSources(); ++src)
+            {
+                const Real* inputExcitation = input
+                                              .getSingleSampleReadPointer
+                                              (src, ear, 0);
+
+                Real* outputExcitation = output_
+                                         .getSingleSampleWritePointer
+                                         (src, ear, 0);
+
+                for (int chn = 0; chn < input.getNChannels(); ++chn)
+                {
+                    Real eNoise = eTot[chn] - inputExcitation[chn];
+                    Real threshold = k_[chn] * eNoise;
+                    if (inputExcitation[chn] > threshold)
+                        outputExcitation[chn] = inputExcitation[chn] - threshold;
+                    else
+                        outputExcitation[chn] = 0.0;
+                }
+            }
+        }
+    }
+
+    //output SignalBanks are cleared so not to worry about filter state
+    void SpecificPartialLoudnessCHGM2011::resetInternal()
+    {
+    }
+}

src/modules/SpecificPartialLoudnessCHGM2011.h

@@ -0,0 +1,45 @@
+/*
+ * Copyright (C) 2014 Dominic Ward <contactdominicward@gmail.com>
+ *
+ * This file is part of Loudness
+ *
+ * Loudness is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Loudness is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Loudness.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+#ifndef SpecificPartialLoudnessCHGM2011_H
+#define SpecificPartialLoudnessCHGM2011_H
+
+#include "../support/Module.h"
+
+namespace loudness{
+
+    class SpecificPartialLoudnessCHGM2011 : public Module
+    {
+    public:
+ 
+        SpecificPartialLoudnessCHGM2011();
+
+        virtual ~SpecificPartialLoudnessCHGM2011();
+
+    private:
+        virtual bool initializeInternal(const SignalBank &input);
+        virtual bool initializeInternal(){return 0;};
+        virtual void processInternal(const SignalBank &input);
+        virtual void processInternal(){};
+        virtual void resetInternal();
+
+        RealVec k_;
+    };
+}
+#endif

src/modules/SpecificPartialLoudnessMGB1997.cpp

@@ -183,7 +183,7 @@ namespace loudness{
                 for (int chn = 0; chn < input.getNChannels(); ++chn)
                 {
                     Real eNoise = eTot[chn] - eSig[chn];
-                    Real eThrn = eNoise * kParam_[chn];
+                    Real eThrn = kParam_[chn] * eNoise + eThrqParam_[chn];
                     Real nSig = 0.0;
 
                     if (eSig[chn] > 1e-10)
@@ -192,29 +192,32 @@ namespace loudness{
                         {
                             if (eSig[chn] >= eThrn) // Equation 19
                             {
-                                nSig = parameterC2_ * std::pow (eTot[chn],
-                                       yearExp_) - parameterC2_
-                                       * (std::pow (eNoise + eThrn 
-                                       + eThrqParam_[chn], yearExp_)
-                                       - std::pow (eThrqParam_[chn] * gParam_[chn]
-                                       + aParam_[chn], alphaParam_[chn])
-                                       + std::pow (aParam_[chn],
-                                       alphaParam_[chn])) * std::pow (eThrn
-                                       / eSig[chn], 0.3);
+                                nSig = parameterC2_ *
+                                       std::pow (eTot[chn], yearExp_);
+
+                                nSig -= parameterC2_ * 
+                                        (std::pow (eNoise + eThrn, yearExp_) -
+                                        std::pow (eThrqParam_[chn] *
+                                        gParam_[chn] + aParam_[chn],
+                                        alphaParam_[chn]) +
+                                        std::pow (aParam_[chn],
+                                        alphaParam_[chn])) *
+                                        std::pow (eThrn / eSig[chn], 0.3);
                             }
                             else // Equation 20
                             {
                                 nSig = parameterC_ * std::pow (2.0 * eSig[chn]
-                                       / (eSig[chn] + eThrn), 1.5)
-                                       * ((std::pow (eThrqParam_[chn]
-                                       * gParam_[chn] + aParam_[chn],
-                                       alphaParam_[chn])
-                                       - std::pow (aParam_[chn], alphaParam_[chn]))
-                                       / (std::pow (eNoise + eThrn
-                                       + eThrqParam_[chn], yearExp_)
-                                       - std::pow(eNoise, yearExp_)))
-                                       * (std::pow (eTot[chn], yearExp_)
-                                       - std::pow (eNoise, yearExp_));
+                                       / (eSig[chn] + eThrn), 1.5);
+
+                                nSig *= (std::pow (eThrqParam_[chn] *
+                                        gParam_[chn] + aParam_[chn],
+                                        alphaParam_[chn]) - 
+                                        std::pow (aParam_[chn],
+                                        alphaParam_[chn])) /
+                                        (std::pow (eNoise + eThrn, yearExp_) -
+                                        std::pow(eNoise, yearExp_));
+                                nSig *= std::pow (eTot[chn], yearExp_)
+                                        - std::pow (eNoise, yearExp_);
                             }
                         }
                         else // eTot <= 100 dB
@@ -224,33 +227,34 @@ namespace loudness{
                                 nSig = parameterC_ * (std::pow (eTot[chn]
                                        * gParam_[chn] + aParam_[chn],
                                        alphaParam_[chn]) - std::pow (
-                                       aParam_[chn], alphaParam_[chn]))
-                                       - parameterC_ * (std::pow (
-                                       gParam_[chn] * (eNoise + eThrn
-                                       + eThrqParam_[chn]) + aParam_[chn],
-                                       alphaParam_[chn]) - std::pow (
-                                       eThrqParam_[chn] * gParam_[chn]
+                                       aParam_[chn], alphaParam_[chn]));
+
+                                nSig -= parameterC_ * (std::pow (
+                                       gParam_[chn] * (eNoise + eThrn) +
+                                       aParam_[chn], alphaParam_[chn]) -
+                                       std::pow (eThrqParam_[chn] * gParam_[chn]
                                        + aParam_[chn], alphaParam_[chn]))
                                        * std::pow (eThrn / eSig[chn], 0.3);
                             }
                             else // Equation 18
                             {
                                 nSig = parameterC_ * std::pow (2.0 * eSig[chn]
-                                       / (eSig[chn] + eThrn), 1.5)
-                                       * ((std::pow (eThrqParam_[chn]
-                                       * gParam_[chn] + aParam_[chn],
-                                       alphaParam_[chn]) - std::pow (
-                                       aParam_[chn], alphaParam_[chn]))
-                                       / (std::pow ((eNoise + eThrn
-                                       + eThrqParam_[chn]) * gParam_[chn]
-                                       + aParam_[chn], alphaParam_[chn])
-                                       - std::pow (eNoise * gParam_[chn]
-                                       + aParam_[chn], alphaParam_[chn])))
-                                       * ( std::pow (eTot[chn]
-                                       * gParam_[chn] + aParam_[chn],
-                                       alphaParam_[chn]) - std::pow (eNoise
-                                       * gParam_[chn] + aParam_[chn],
-                                       alphaParam_[chn]));
+                                       / (eSig[chn] + eThrn), 1.5);
+                                
+                                nSig *= (std::pow (eThrqParam_[chn] *
+                                        gParam_[chn] + aParam_[chn],
+                                        alphaParam_[chn]) - std::pow (
+                                        aParam_[chn], alphaParam_[chn])) /
+                                        (std::pow ((eNoise + eThrn) *
+                                        gParam_[chn] + aParam_[chn],
+                                        alphaParam_[chn]) - std::pow (
+                                        eNoise * gParam_[chn] + aParam_[chn],
+                                        alphaParam_[chn]));
+
+                                nSig *= std::pow (eTot[chn] * gParam_[chn] +
+                                        aParam_[chn], alphaParam_[chn]) -
+                                        std::pow (eNoise * gParam_[chn] +
+                                        aParam_[chn], alphaParam_[chn]);
                             }
                         }
                     }

src/python/__init__.py

@@ -1,5 +1,5 @@
-import sound
-import extractors
-import iterators
-import predictors
-import spectra
+from . import sound
+from . import extractors
+from . import iterators
+from . import predictors
+from . import spectra

src/python/extractors.py

@@ -7,7 +7,7 @@
 
 class StationaryLoudnessExtractor:
 
-    def __init__(self, model, outputs):
+    def __init__(self, model, outputs, alwaysReinitialize=True):
 
         if model.isDynamic():
             raise ValueError("Model cannot be dynamic")
@@ -23,6 +23,7 @@ def __init__(self, model, outputs):
         self.nEars = 0
         self.bank = ln.SignalBank()
         self.initialize = True
+        self.alwaysReinitialize = alwaysReinitialize
 
     def process(self, frequencies, intensityLevels):
 
@@ -51,7 +52,8 @@ def process(self, frequencies, intensityLevels):
             self.outputDict['Frequencies'] = frequencies
             self.nEars = nEars
             self.model.initialize(self.bank)
-        self.initialize = False
+        if not self.alwaysReinitialize:
+            self.initialize = False
 
         self.bank.setSignals(intensities.T.reshape((1, nEars, nComponents, 1)))
 
@@ -121,120 +123,176 @@ def __init__(self,
                                  self.hopSize,
                                  self.fs)
         self.outputs = outputs
-        self.gainInDecibels = gainInDecibels
         if self.outputs is not None:
-            if type(self.outputs) is not list:
+            if isinstance(self.outputs, str):
                 self.outputs = [self.outputs]
-            self.model.setOutputsToAggregate(self.outputs)
         else:
             raise ValueError("Must specify outputs")
 
+        '''
+        if aggregate:
+            self.model.setOutputsToAggregate(self.outputs)
+        '''
+
         if not self.model.initialize(self.inputBuf):
             raise ValueError("Problem initialising the model!")
 
         self.outputDict = {}
-        self.nSamplesToPadStart = np.round(numSecondsToPadStartBy * self.fs)
-        self.nSamplesToPadEnd = np.round(numSecondsToPadEndBy * self.fs)
+        self.nSamplesToPadStart = int(
+                np.round(numSecondsToPadStartBy * self.fs))
+        self.nSamplesToPadEnd = int(
+                np.round(numSecondsToPadEndBy * self.fs))
         self.frameTimeOffset = frameTimeOffset
         self.x = None
         self.processed = False
         self.loudness = None
         self.globalLoudness = None
 
-    def process(self, inputSignal):
-        '''
-        Process the numpy array `inputSignal' using a dynamic loudness
-        model.  For stereo signals, the input signal must have two dimensions
-        with shape (nSamples x 2).  For monophonic signals, the input signal
-        can be 2D, i.e. (nSamples x 1), or one dimensional. For multiple
-        sources, enter a list of inputSignals (easier than faffing with
-        reshaping).
-        '''
-        # Input checks
-        if type(inputSignal) is not list:
+    def configureInput(self, inputSignal):
+
+        if not isinstance(inputSignal, list):
             if inputSignal.ndim == 1:
-                self.inputSignal = inputSignal.reshape(1,
-                                                       1,
-                                                       1,
-                                                       inputSignal.size)
+                sig = inputSignal.reshape(1,
+                                          1,
+                                          1,
+                                          inputSignal.size)
             elif inputSignal.ndim < 3:
-                self.inputSignal = inputSignal.T.reshape(1,
-                                                         inputSignal.shape[1],
-                                                         1,
-                                                         inputSignal.shape[0])
+                sig = inputSignal.T.reshape(1,
+                                            inputSignal.shape[1],
+                                            1,
+                                            inputSignal.shape[0])
         else:
             maxSamples = np.max([x.shape[0] for x in inputSignal])
-            self.inputSignal = np.zeros((self.nInputSources,
-                                         self.nInputEars,
-                                         1,
-                                         maxSamples))
-            for i, sig in enumerate(inputSignal):
-                self.inputSignal[i, :, 0, 0:sig.shape[0]] = sig.T
-
-        '''Pad end so that we can obtain analysis over the last sample,
-        No neat way to do this at the moment so assume 0.2ms is enough.'''
-        self.inputSignal = np.concatenate((
+            sig = np.zeros((self.nInputSources,
+                            self.nInputEars,
+                            1,
+                            maxSamples))
+            for i, isig in enumerate(inputSignal):
+                sig[i, :, 0, :isig.shape[0]] = isig.T
+
+        sig = np.concatenate((
             np.zeros((self.nInputSources, self.nInputEars,
                       1, self.nSamplesToPadStart)),
-            self.inputSignal,
+            sig,
             np.zeros((self.nInputSources, self.nInputEars,
                       1, self.nSamplesToPadEnd))), 3)
 
-        # Apply any gain
-        self.inputSignal *= 10 ** (self.gainInDecibels / 20.0)
-
         # configure the number of output frames needed
         nOutputFrames = int(
-            np.ceil(self.inputSignal.shape[3] / float(self.hopSize))
+            np.ceil(sig.shape[3] / float(self.hopSize))
         )
-        self.outputDict['FrameTimes'] = (
-            self.frameTimeOffset + np.arange(nOutputFrames) *
-            self.hopSize / float(self.fs)
+        return sig, nOutputFrames
+
+    def outputToHDF5(self, inputSignal, nOutputFrames, hdf5Group):
+
+        hdf5Group.clear()
+        hdf5Group.create_dataset(
+            'FrameTime',
+            data = self.frameTimeOffset + np.arange(nOutputFrames) *
+            self.hopSize / float(self.fs),
         )
 
+        outputBanks = []
+        datasets = []
+        for name in self.outputs:
+            bank = self.model.getOutput(name)
+            outputBanks.append(bank)
+            shape = (
+                nOutputFrames,
+                bank.getNSources(),
+                bank.getNEars(),
+                bank.getNChannels(),
+                bank.getNSamples(),
+            )
+
+            datasets.append(hdf5Group.create_dataset(name, shape, chunks=True))
+
         # One process call every hop samples
         for frame in range(nOutputFrames):
 
-            # Any overlap is generated c++ side, so just fill the input
-            # SignalBank with new blocks
             startIdx = frame * self.hopSize
             endIdx = startIdx + self.hopSize
             self.inputBuf.setSignals(
-                self.inputSignal[:, :, :, startIdx:endIdx]
+                inputSignal[:, :, :, startIdx:endIdx]
             )
 
             # Process the input buffer
             self.model.process(self.inputBuf)
 
-        # get outputs
+            # Store
+            for bank, dataset in zip(outputBanks, datasets):
+                dataset[frame] = bank.getSignals()
+
+        # Processing complete so clear internal states
+        self.model.reset()
+        self.processed = True
+
+    def outputToDictionary(self, inputSignal, nOutputFrames):
+
+        dic = {'FrameTime':
+               self.frameTimeOffset + np.arange(nOutputFrames) *
+               self.hopSize / float(self.fs),
+        }
+
+        outputBanks = []
         for name in self.outputs:
             bank = self.model.getOutput(name)
-            self.outputDict[name] = np.copy(
-                np.squeeze(bank.getAggregatedSignals())
+            outputBanks.append(bank)
+            shape = (
+                nOutputFrames,
+                bank.getNSources(),
+                bank.getNEars(),
+                bank.getNChannels(),
+                bank.getNSamples(),
             )
-            if self.outputDict[name].ndim == 1:
-                self.outputDict[name] = self.outputDict[name].reshape((-1, 1))
+            dic[name] = np.squeeze(np.zeros(shape=shape))
+
+        # One process call every hop samples
+        for frame in range(nOutputFrames):
+
+            startIdx = frame * self.hopSize
+            endIdx = startIdx + self.hopSize
+            self.inputBuf.setSignals(
+                inputSignal[:, :, :, startIdx:endIdx]
+            )
+
+            # Process the input buffer
+            self.model.process(self.inputBuf)
+
+            # Store
+            for bank, name in zip(outputBanks, self.outputs):
+                dic[name][frame] = np.squeeze(bank.getSignals())
+
+        return dic
+
+
+    def process(self, inputSignal, hdf5Group=None):
+        '''
+        Process the numpy array `inputSignal' using a dynamic loudness
+        model.  For stereo signals, the input signal must have two dimensions
+        with shape (nSamples x 2).  For monophonic signals, the input signal
+        can be 2D, i.e. (nSamples x 1), or one dimensional. For multiple
+        sources, enter a list of inputSignals (easier than faffing with
+        reshaping).
+        '''
+
+        if hdf5Group is None:
+            sig, nOutputFrames = self.configureInput(inputSignal)
+            dic = self.outputToDictionary(sig, nOutputFrames)
+        elif isinstance(hdf5Group, h5py.Group):
+            sig, nOutputFrames = self.configureInput(inputSignal)
+            self.outputToHDF5(sig, nOutputFrames, hdf5Group)
 
         # Processing complete so clear internal states
         self.model.reset()
         self.processed = True
 
+        if hdf5Group is None:
+            return dic
+
     def reinitializeModel(self):
         self.model.initialize(self.inputBuf)
 
-    def saveOutputToPickle(self, filename):
-        '''
-        Saves the complete output dictionary to a pickle file.
-        '''
-        if self.processed:
-            print 'Saving to pickle file'
-            filename, ext = os.path.splitext(filename)
-            with open(filename + '.pickle', 'wb') as outfile:
-                pickle.dump(
-                    self.outputDict,
-                    outfile,
-                    protocol=pickle.HIGHEST_PROTOCOL)
-
 
 class BatchWavFileProcessor:
     """Class for processing multiple wav files using a given loudness model.
@@ -303,7 +361,7 @@ def process(self, model):
 
         model.setOutputsToAggregate(self.outputs)
 
-        print "Output will be saved to ", self.filename
+        print ("Output will be saved to ", self.filename)
         h5File = h5py.File(self.filename, 'w')
 
         processor = ln.AudioFileProcessor(
@@ -315,7 +373,7 @@ def process(self, model):
 
         for wavFile in self.wavFiles:
 
-            print("Processing file %s ..." % wavFile)
+            print ("Processing file %s ..." % wavFile)
 
             # Create new group
             wavFileGroup = h5File.create_group(wavFile)

src/python/iterators.py

@@ -1,8 +1,8 @@
 import os
 import numpy as np
 import matplotlib.pyplot as plt
-from sound import Sound
-from extractors import DynamicLoudnessExtractor, StationaryLoudnessExtractor
+from .sound import Sound
+from .extractors import DynamicLoudnessExtractor, StationaryLoudnessExtractor
 
 
 def asIs(x):
@@ -33,7 +33,8 @@ def __init__(self,
                  nIters=10,
                  alpha=1.0,
                  nSecondsToPadStartBy=0.0,
-                 nSecondsToPadEndBy=0.2):
+                 nSecondsToPadEndBy=0.2,
+                 printResults=False):
 
         if type(outputName) is list:
             outputName = outputName[0]
@@ -50,7 +51,7 @@ def __init__(self,
         self.tol = tol
         self.nIters = nIters
         self.alpha = alpha
-        self.printResults = False
+        self.printResults = printResults
 
         if loudnessFunction is None:
             self.loudnessFunction = np.mean
@@ -59,9 +60,9 @@ def __init__(self,
 
     def extractLoudness(self, signal, gainInDecibels=0.0):
 
-        self.extractor.gainInDecibels = gainInDecibels
-        self.extractor.process(signal)
-        timeSeries = self.extractor.outputDict[self.outputName]
+        gain = 10 ** (gainInDecibels / 20.0)
+        out = self.extractor.process(signal * gain)
+        timeSeries = out[self.outputName]
         loudness = self.loudnessFunction(timeSeries)
         return loudness
 
@@ -107,10 +108,12 @@ def __init__(self,
                  loudnessFunction=None,
                  tol=0.1,
                  nIters=10,
-                 alpha=1.0):
+                 alpha=1.0,
+                 alwaysReinitialize=False):
 
         self.outputName = outputName
-        self.extractor = StationaryLoudnessExtractor(model, outputName)
+        self.extractor = StationaryLoudnessExtractor(model, outputName,
+                                                     alwaysReinitialize)
         self.converged = False
         self.tol = tol
         self.nIters = nIters
@@ -184,29 +187,34 @@ def __init__(self,
                  nIters=10,
                  alpha=1.0,
                  nSecondsToPadStartBy=0.0,
-                 nSecondsToPadEndBy=0.2):
+                 nSecondsToPadEndBy=0.2,
+                 printResults=False):
 
         self.nSecondsToPadStartBy = nSecondsToPadStartBy
         self.nSecondsToPadEndBy = nSecondsToPadEndBy
         self.loudnessFunction = loudnessFunction
         self.tol = tol
         self.nIters = nIters
         self.alpha = alpha
+        self.model = model
+        self.fs = fs
+        self.printResults = printResults
 
         self.output = output
         if self.output is None or type(self.output) is list:
             raise ValueError("Only one model output allowed.")
 
-    def process(self, listOfInputSignals, targetLoudness):
+    def process(self, listOfInputSignals, targetLoudness=None):
+        print(self.printResults)
 
         hasMultipleTargets = False
-        if type(self.targetLoudness) in [np.ndarray, list]:
-            if len(self.targetLoudness) == len(listOfInputSignals):
+        if type(targetLoudness) in [np.ndarray, list]:
+            if len(targetLoudness) == len(listOfInputSignals):
                 hasMultipleTargets = True
 
-        self.gains = np.zeros(len(self.listOfInputSignals))
+        gains = np.zeros(len(listOfInputSignals))
 
-        for i, signal in enumerate(self.listOfInputSignals):
+        for i, signal in enumerate(listOfInputSignals):
 
             if signal.ndim == 1:
                 signal = signal.reshape((-1, 1))
@@ -222,22 +230,23 @@ def process(self, listOfInputSignals, targetLoudness):
                     self.nIters,
                     self.alpha,
                     self.nSecondsToPadStartBy,
-                    self.nSecondsToPadEndBy)
+                    self.nSecondsToPadEndBy,
+                    self.printResults)
 
             # Process the audio file
             print("Processing signal %d ..." % i)
 
             if hasMultipleTargets:
-                target = self.targetLoudness[i]
+                target = targetLoudness[i]
             else:
-                target = self.targetLoudness
+                target = targetLoudness
 
             gain = processor.process(signal, target)
 
             # Get gain required for target loudness
-            self.gains[i] = gain
+            gains[i] = gain
 
-        return self.gains
+        return gains
 
 
 class EqualDynamicLoudnessIterator():

src/python/predictors.py

@@ -1,9 +1,9 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from sound import Sound
-from iterators import StationaryLoudnessIterator, DynamicLoudnessIterator
-from extractors import StationaryLoudnessExtractor, DynamicLoudnessExtractor
-from scipy.interpolate import interp1d
+from .sound import Sound
+from .iterators import StationaryLoudnessIterator, DynamicLoudnessIterator
+from .extractors import StationaryLoudnessExtractor, DynamicLoudnessExtractor
+from  scipy.interpolate import interp1d
 
 # ISO 389-7 - free-field values
 freqsISO389 = np.array([20.0, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200,
@@ -101,7 +101,7 @@ def computeErrors(self):
             maxE = np.max(np.abs(error))
             return (error, rMSE, maxE)
         else:
-            print 'Iterative process did not converge.'
+            print ('Iterative process did not converge.')
             return (None, None, None)
 
     def getResults(self):
@@ -219,7 +219,7 @@ def computeErrors(self):
             maxE = np.max(np.abs(error))
             return (error, rMSE, maxE)
         else:
-            print 'Iterative process did not converge.'
+            print ('Iterative process did not converge.')
             return (None, None, None)
 
     def getResults(self):

src/python/sound.py

@@ -188,7 +188,7 @@ def applyRamp(self, dur):
         theta = np.pi * 0.5 * np.arange(d) / float(d)
         nSamplesmiddle = self.nSamples - 2 * d
         if nSamplesmiddle < 0:
-            print "Duration too large"
+            print ("Duration too large")
         else:
             self.data *= np.concatenate((
                 np.cos(theta + np.pi * 1.5),

src/support/AuditoryTools.cpp

@@ -163,10 +163,10 @@ namespace loudness{
     void OME::getData()
     {
         //The data
-        Real freqs[41] = {0, 20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250,
-            315, 400, 500, 630, 750, 800, 1000, 1250, 1500, 1600, 2000, 2500,
-            3000, 3150, 4000, 5000, 6000, 6300, 8000, 9000, 10000, 11200, 12500,
-            14000, 15000, 16000, 18000, 20000};
+        Real freqs[41] = {0, 20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200,
+            250, 315, 400, 500, 630, 750, 800, 1000, 1250, 1500, 1600, 2000,
+            2500, 3000, 3150, 4000, 5000, 6000, 6300, 8000, 9000, 10000, 11200,
+            12500, 14000, 15000, 16000, 18000, 20000};
 
         Real freeANSI[41] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0.1, 0.3, 0.5, 0.9, 1.4,
                 1.6, 1.7, 2.5, 2.7, 2.6, 2.6, 3.2, 5.2, 6.6, 12, 16.8, 15.3, 15.2,
@@ -176,11 +176,11 @@ namespace loudness{
             1.6, 1.7, 2.2, 2.7, 2.9, 3.8, 5.3, 6.8, 7.2, 10.2, 14.9, 14.5, 14.4, 
             12.7, 10.8, 8.9, 8.7, 8.5, 6.2, 5, 4.5, 4, 3.3, 2.6, 2, 2, 2};
 
-        Real midANSI[41] = {-100, -39.6, -32, -25.85, -21.4, -18.5, -15.9, -14.1,
-            -12.4, -11, -9.6, -8.3, -7.4, -6.2, -4.8, -3.8, -3.3, -2.9, 
-            -2.6, -2.6, -4.5, -5.4, -6.1, -8.5,-10.4, -7.3, -7, -6.6,
-            -7, -9.2, -10.2, -12.2, -10.8, -10.1, -12.7, -15, -18.2,
-            -23.8, -32.3, -45.5, -50};
+        Real midANSI[41] = {-100, -39.6, -32, -25.85, -21.4, -18.5, -15.9,
+            -14.1, -12.4, -11, -9.6, -8.3, -7.4, -6.2, -4.8, -3.8, -3.3, -2.9,
+            -2.6, -2.6, -4.5, -5.4, -6.1, -8.5,-10.4, -7.3, -7, -6.6, -7, -9.2,
+            -10.2, -12.2, -10.8, -10.1, -12.7, -15, -18.2, -23.8, -32.3, -45.5,
+            -50};
 
         Real midWard[41] = {-36.33, -36.33, -31.7, -27.24, -23.08, -19.61,
             -16.41, -13.4, -10.88, -8.67, -6.6, -4.88, -3.46, -2.23, -0.41,
@@ -196,25 +196,39 @@ namespace loudness{
             -21.8, -22.3, -23.2, -24};
 
         //Beyerdynamic DT990 requires own sample points
-        Real phoneFreqs[86] = {0, 5, 11, 17, 23, 29, 35, 41, 46, 52, 58, 64, 70,
-            76, 82, 88, 101, 116, 133, 153, 176, 202, 232, 266, 306, 352, 404,
-            464, 533, 612, 703, 808, 928, 1066, 1224, 1406, 1615, 1856, 2131,
-            2448, 2813, 3231, 3711, 3928, 4144, 4396, 4897, 5215, 5437, 5625,
-            5780, 6110, 6500, 6800, 7422, 7920, 8400, 8526, 8830, 9000, 9160,
-            9230, 9360, 9650, 9900, 10500, 10900, 11250, 11700, 12200, 12500,
-            12923, 13500, 14000, 14350, 14600, 14750, 14845, 15350, 16100,
-            16700, 17200, 17750, 18420, 19587, 20000};
-
-        Real dt990[86] = {-6.45, -6.45, -7.48, -5.75, -3.69, -2.52, -1.61,
-            -0.71, -0.42, 0.13,  0.24, -0.51,  0.56,  1.61,  2.02,  2.28,  2.6 ,
-            2.79, 2.82,  2.69,  2.43,  2.09,  1.68,  1.31,  0.77,  0.15, -0.31,
-            -0.64, -0.96, -1.36, -1.41, -1.16, -0.52,  0.  ,  1.24,  1.84, 2.8 ,
-            4.28,  5.77,  7.5 ,  8.39,  7.2 ,  4.3 ,  3.37,  4.  , 4.52,  2.25,
-            1.36,  2.26,  4.69,  7.01,  9.37,  9.6 ,  9.66, 7.8 ,  6.07,  8.25,
-            8.9 ,  7.84,  5.42,  2.24,  0.73, -0.08, 0.69,  1.23,  2.01,  1.47,
-            0.8 ,  0.54,  3.06,  3.5 ,  1.79, -3.52, -7.86, -8.68, -9.81, -7.1 ,
-            -5.09, -1.14,  0.02,  1.6 , 0.25, -2.8 , -3.97, -7.45, -6.8 };
-            
+        Real phoneFreqs[88] = {0, 5, 11, 17, 23,    29,    35,
+            41,    46, 52,    58,    64,    70,    76,    82,    88,   101,
+            116, 133,   153,   176,   202,   232,   266,   306,   352,   404,
+            464,   533,   612,   703,   808,   928,  1000,  1224,  1406, 1615,
+            1856,  2131,  2448,  2813,  3231,  3711,  3928,  4144, 4396, 4897,
+            5215,  5437,  5625,  5780,  5900, 6100, 6500, 7000, 7422,  7920,
+            8400, 8526, 8830, 9000,  9160,9230, 9360, 9650, 9900, 10500, 10900,
+            11300, 11700, 12000, 12200, 12700,13500,14000, 14350, 14600,14750,
+            14845, 15350,16100, 16700,17000, 17200, 17750, 18420,19587,20000};
+
+        Real dt990[88] = {-6.40, -6.40, -7.62, -5.55, -3.56, -2.47, -1.60,
+            -0.72, -0.39, 0.20, 0.22, -0.57, 0.61, 1.62, 2.00, 2.25, 2.56,
+            2.78, 2.80, 2.69, 2.43, 2.13, 1.63, 1.32, 0.78, 0.14, -0.33, -0.64,
+            -0.98, -1.37, -1.43, -1.18, -0.56, 0.00, 1.24, 1.84, 2.82, 4.31,
+            5.82, 7.50, 8.37, 7.20, 4.32, 3.49, 3.99, 4.52, 2.26, 1.36, 2.25,
+            4.67, 6.95, 8.37, 9.34, 9.63, 9.53, 7.98, 6.36, 8.78, 9.49, 8.85,
+            6.96, 4.83, 3.87, 2.43, 0.83, 1.47, 2.81, 1.55, 0.74, 0.73, 1.97,
+            3.25, 3.10, -3.52, -7.35, -8.66, -5.70, -3.47, -2.35, -0.35, 0.09,
+            1.58, 1.29, 0.28, -2.79, -3.99, -7.43, -6.36};
+
+        // Based on MIT KEMAR ELEVATION=0 AZIMUTH=0
+        Real kemarAvg00Freqs[42] = {100, 260, 300, 452, 581, 775, 904, 1080,
+            1270, 1572, 1852, 2239, 2476, 2713, 2900, 3252, 3600, 4100, 4300,
+            5017, 5211, 5491, 5900, 6072, 6309, 7321, 8300, 9173, 9539, 9755,
+            10000, 10700, 11370, 12037, 13200, 13717, 15159, 16215, 16753,
+            17162, 19251, 20000};
+
+        Real kemarAvg00[42] = {0.00, 1.00, 1.13, 0.99, 2.00, 2.02, 1.56, 2.22,
+            2.03, 5.35, 10.24, 13.77, 12.80, 10.73, 10.08, 9.59, 10.30, 10.35,
+            8.49, 5.98, 6.49, 6.57, 4.27, 2.75, 1.05, -3.83, -10.63, -5.24,
+            -2.34, 0.34, 3.54, 4.24, 2.54, 2.92, 9.30, 10.01, 4.80, -4.14,
+            -3.76, -4.44, -7.55, -12.19};
+
         //frequencies to use for middle ear interpolation
         middleEarFreqPoints_.assign(freqs, freqs + 41);
 
@@ -249,8 +263,13 @@ namespace loudness{
                 outerEarFreqPoints_ = middleEarFreqPoints_;
                 break;
             case BD_DT990:
-                outerEardB_.assign(dt990, dt990 + 86);
-                outerEarFreqPoints_.assign(phoneFreqs, phoneFreqs + 86);
+                outerEardB_.assign(dt990, dt990 + 88);
+                outerEarFreqPoints_.assign(phoneFreqs, phoneFreqs + 88);
+                break;
+            case LML_FREEFIELD:
+                outerEardB_.assign(kemarAvg00, kemarAvg00 + 42);
+                outerEarFreqPoints_.assign(kemarAvg00Freqs,
+                                           kemarAvg00Freqs + 42);
                 break;
             default:
                 outerEarType_ = NONE;
@@ -307,13 +326,17 @@ namespace loudness{
             //outer ear
             if(outerEarType_ != NONE)
             {
+                Real firstFreq = outerEarFreqPoints_.front();
+                Real firstDataPoint = outerEardB_.front();
                 Real lastFreq = outerEarFreqPoints_.back();
                 Real lastDataPoint = outerEardB_.back();
                 s.set_points(outerEarFreqPoints_, outerEardB_);
 
                 for(uint i = 0; i < freqs.size(); i++)
                 {
-                    if(freqs[i] >= lastFreq)
+                    if(freqs[i] <= firstFreq)
+                        response_[i] += firstDataPoint;
+                    else if(freqs[i] >= lastFreq)
                         response_[i] += lastDataPoint;
                     else
                         response_[i] += s(freqs[i]);

src/support/AuditoryTools.h

@@ -65,6 +65,16 @@ namespace loudness{
         return (pow (10.0, (cam / 21.366)) - 1.0) / 4368e-6;
     }
 
+    inline Real hertzToNormalisedCam (Real freq)
+    {
+        return (hertzToCam(freq) / hertzToCam(1000.0)) - 1;
+    }
+
+    inline Real normalisedCamToHertz (Real normalisedCam)
+    {
+        return (normalisedCam + 1) * hertzToCam(1000.0);
+    }
+
     /**
      * @brief Returns the ERB in Bark of the auditory filter for normally hearing
      * listeners at @a centreFreq in Hz.
@@ -213,7 +223,8 @@ namespace loudness{
                 WARD_MIDDLE_EAR,
                 ANSIS342007_FREEFIELD,
                 ANSIS342007_DIFFUSEFIELD,
-                BD_DT990
+                BD_DT990,
+                LML_FREEFIELD
             };
 
             OME(){};

src/support/Common.h

@@ -22,11 +22,11 @@
 
 #include <memory>
 #include <vector>
+#include <cmath>
 #include <string>
 #include <algorithm>
 #include <map>
 #include "Debug.h"
-#include "UsefulFunctions.h"
 
 /*
  * Expose objects

src/support/Filter.cpp

@@ -22,12 +22,10 @@
 
 namespace loudness{
 
-    Filter::Filter() : gain_(1.0)
-    {}
+    Filter::Filter() {}
 
     //Is this needed?
     Filter::Filter(int order) :
-        gain_(1.0),
         order_(order)
     {}
 
@@ -118,18 +116,9 @@ namespace loudness{
         return aCoefs_;
     }
 
-    Real Filter::getGain() const
-    {
-        return gain_;
-    }
-
     int Filter::getOrder() const
     {
         return order_;
     }
 
-    void Filter::setGain(Real gain)
-    {
-        gain_ = gain;
-    }
 }

src/support/Filter.h

@@ -79,16 +79,6 @@ namespace loudness{
         const RealVec& getBCoefs() const;
         const RealVec& getACoefs() const;
 
-        /**
-         * @brief Sets the linear gain of the digital filter.
-         */
-        void setGain(Real gain);
-
-        /**
-         * @brief Returns the linear gain of the digital filter.
-         */
-        Real getGain() const;
-
         /**
          * @brief Returns the order of the digital filter.
          *

src/support/SignalBank.h

@@ -21,6 +21,7 @@
 #define SIGNALBANK_H
 
 #include "Common.h"
+#include "UsefulFunctions.h"
 
 namespace loudness{
 
@@ -243,7 +244,7 @@ namespace loudness{
         /** Get the value of a single sample in a given ear and channel. 
          * watch your bounds.
          */
-        inline Real getSample(int source, int ear, int channel, int sample) const
+        inline Real getSample(int source=0, int ear=0, int channel=0, int sample=0) const
         {
             LOUDNESS_ASSERT(
                     isPositiveAndLessThanUpper(source, nSources_) &&

src/support/UsefulFunctions.h

@@ -20,13 +20,19 @@
 #ifndef  USEFULFUNCTIONS_H
 #define  USEFULFUNCTIONS_H
 
-#include <cmath>
-#include <vector>
+#include "Common.h"
 #define PI 3.14159265358979323846264338327
 #define LOWER_LIMIT_DB -100
 
 namespace loudness{
 
+    inline void killDenormal (Real& value)
+    {
+        static const Real antiDenormal = 1e-18;
+        value += antiDenormal;
+        value -= antiDenormal;
+    }
+
     /** Factorial of a number. */
     template <typename Type>
     inline Type factorial (Type const& value)