ENH: noise models - added A, B, C, and D weightings to SPL

The SPL can now be weighted according to the new 'SPLweighting' entry
that can be set to:
- none: no weighting
- dBA : dB(A)
- dBB : dB(B)
- dBC : dB(C)
- dBD : dB(D)

This commit also includes code refactoring of the  noiseModel class to
remove the dependency on noiseFFT/declutter.

applications/utilities/postProcessing/noise/Make/options

@@ -1,7 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/randomProcesses/lnInclude
+    -I$(LIB_SRC)/randomProcesses/lnInclude \
+    -I$(FFTW_INC_DIR)
 
 EXE_LIBS = \
     -lsampling \

applications/utilities/postProcessing/noise/noise.C

@@ -87,7 +87,6 @@ Usage
 
 
 See also
-    noiseFFT.H
     noiseModel.H
     windowModel.H
 

src/randomProcesses/noise/noiseFFT/noiseFFT.C

@@ -127,7 +127,7 @@ void Foam::noiseFFT::octaveBandInfo
         }
     }
 
-    fBandIDs = bandIDs.toc();
+    fBandIDs = bandIDs.sortedToc();
 
     if (fc.size())
     {

src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2018 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,8 @@ License
 #include "noiseModel.H"
 #include "functionObject.H"
 #include "argList.H"
+#include "fft.H"
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -37,6 +39,86 @@ namespace Foam
     defineRunTimeSelectionTable(noiseModel, dictionary);
 }
 
+
+const Foam::Enum<Foam::noiseModel::weightingType>
+Foam::noiseModel::weightingTypeNames_
+({
+    {weightingType::none, "dB"},
+    {weightingType::dBA, "dBA"},
+    {weightingType::dBB, "dBB"},
+    {weightingType::dBC, "dBC"},
+    {weightingType::dBD, "dBD"},
+});
+
+
+void Foam::noiseModel::setOctaveBands
+(
+    const scalarField& f,
+    const scalar fLower,
+    const scalar fUpper,
+    const scalar octave,
+    labelList& fBandIDs,
+    scalarField& fCentre
+)
+{
+    // Set the indices of to the lower frequency bands for the input frequency
+    // range. Ensure that the centre frequency passes though 1000 Hz
+
+    // Low frequency bound given by:
+    //     fLow = f0*(2^(0.5*bandI/octave))
+
+    // Initial (lowest centre frequency)
+    scalar fTest = 15.625;
+
+    const scalar fRatio = pow(2, 1.0/octave);
+    const scalar fRatioL2C = pow(2, 0.5/octave);
+
+    // IDs of band IDs
+    labelHashSet bandIDs(f.size());
+
+    // Centre frequencies
+    DynamicList<scalar> fc;
+
+    // Convert to lower band limit
+    fTest /= fRatioL2C;
+
+    forAll(f, i)
+    {
+        if (f[i] >= fTest)
+        {
+            // Advance band if appropriate
+            while (f[i] > fTest)
+            {
+                fTest *= fRatio;
+            }
+            fTest /= fRatio;
+
+            if (bandIDs.insert(i))
+            {
+                // Also store (next) centre frequency
+                fc.append(fTest*fRatioL2C);
+            }
+            fTest *= fRatio;
+
+            if (fTest > fUpper)
+            {
+                break;
+            }
+        }
+    }
+
+    fBandIDs = bandIDs.sortedToc();
+
+    if (fc.size())
+    {
+        // Remove the last centre frequency (beyond upper frequency limit)
+        fc.remove();
+
+        fCentre.transfer(fc);
+    }
+}
+
+
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 void Foam::noiseModel::readWriteOption
@@ -71,15 +153,13 @@ Foam::scalar Foam::noiseModel::checkUniformTimeStep
 
     if (times.size() > 1)
     {
-        for (label i = 1; i < times.size(); i++)
+        // Uniform time step
+        deltaT = (times.last() - times.first())/scalar(times.size() - 1);
+
+        for (label i = 1; i < times.size(); ++i)
         {
             scalar dT = times[i] - times[i-1];
 
-            if (deltaT < 0)
-            {
-                deltaT = dT;
-            }
-
             if (mag(dT/deltaT - 1) > 1e-8)
             {
                 FatalErrorInFunction
@@ -155,6 +235,267 @@ Foam::fileName Foam::noiseModel::baseFileDir(const label dataseti) const
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::noiseModel::uniformFrequencies
+(
+    const scalar deltaT
+) const
+{
+    const auto& window = windowModelPtr_();
+    const label N = window.nSamples();
+
+    auto tf = tmp<scalarField>::New(N/2 + 1, Zero);
+    auto& f = tf.ref();
+
+    const scalar deltaf = 1.0/(N*deltaT);
+    forAll(f, i)
+    {
+        f[i] = i*deltaf;
+    }
+
+    return tf;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::octaves
+(
+    const scalarField& data,
+    const scalarField& f,
+    const labelUList& freqBandIDs
+) const
+{
+    auto toctData = tmp<scalarField>::New(freqBandIDs.size() - 1, Zero);
+
+    if (freqBandIDs.size() < 2)
+    {
+        WarningInFunction
+            << "Octave frequency bands are not defined "
+            << "- skipping octaves calculation"
+            << endl;
+
+        return toctData;
+    }
+
+    auto& octData = toctData.ref();
+
+    for (label bandI = 0; bandI < freqBandIDs.size() - 1; ++bandI)
+    {
+        label fb0 = freqBandIDs[bandI];
+        label fb1 = freqBandIDs[bandI+1];
+
+        if (fb0 == fb1) continue;
+
+        for (label freqI = fb0; freqI < fb1; ++freqI)
+        {
+            label f0 = f[freqI];
+            label f1 = f[freqI + 1];
+            scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
+            octData[bandI] += dataAve*(f1 - f0);
+        }
+    }
+
+    return toctData;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::Pf(const scalarField& p) const
+{
+    if (planInfo_.active)
+    {
+        if (p.size() != planInfo_.windowSize)
+        {
+            FatalErrorInFunction
+                << "Expected pressure data to have " << planInfo_.windowSize
+                << " values, but received " << p.size() << " values"
+                << abort(FatalError);
+        }
+
+        List<double>& in = planInfo_.in;
+        const List<double>& out = planInfo_.out;
+        forAll(in, i)
+        {
+            in[i] = p[i];
+        }
+
+        ::fftw_execute(planInfo_.plan);
+
+        const label n = planInfo_.windowSize;
+        const label nBy2 = n/2;
+        auto tresult = tmp<scalarField>::New(nBy2 + 1);
+        auto& result = tresult.ref();
+
+        // 0 th value = DC component
+        // nBy2 th value = real only if n is even
+        result[0] = out[0];
+        for (label i = 1; i <= nBy2; ++i)
+        {
+            const auto re = out[i];
+            const auto im = out[n - i];
+            result[i] = sqrt(re*re + im*im);
+        }
+
+        return tresult;
+    }
+
+    return mag(fft::realTransform1D(p));
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::meanPf(const scalarField& p) const
+{
+    const auto& window = windowModelPtr_();
+    const label N = window.nSamples();
+    const label nWindow = window.nWindow();
+
+    auto tmeanPf = tmp<scalarField>::New(N/2 + 1, Zero);
+    auto& meanPf = tmeanPf.ref();
+
+    for (label windowI = 0; windowI < nWindow; ++windowI)
+    {
+        meanPf += Pf(window.apply<scalar>(p, windowI));
+    }
+
+    meanPf /= scalar(nWindow);
+
+    return tmeanPf;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::RMSmeanPf
+(
+    const scalarField& p
+) const
+{
+    const auto& window = windowModelPtr_();
+    const label N = window.nSamples();
+    const label nWindow = window.nWindow();
+
+    scalarField RMSMeanPf(N/2 + 1, Zero);
+    for (label windowI = 0; windowI < nWindow; ++windowI)
+    {
+        RMSMeanPf += sqr(Pf(window.apply<scalar>(p, windowI)));
+    }
+
+    return sqrt(RMSMeanPf/scalar(nWindow))/scalar(N);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::PSDf
+(
+    const scalarField& p,
+    const scalar deltaT
+) const
+{
+    const auto& window = windowModelPtr_();
+    const label N = window.nSamples();
+    const label nWindow = window.nWindow();
+
+    auto tpsd = tmp<scalarField>::New(N/2 + 1, Zero);
+    auto& psd = tpsd.ref();
+
+    for (label windowI = 0; windowI < nWindow; ++windowI)
+    {
+        psd += sqr(Pf(window.apply<scalar>(p, windowI)));
+    }
+
+    scalar fs =  1.0/deltaT;
+    psd /= scalar(nWindow)*fs*N;
+
+    // Scaling due to use of 1-sided FFT
+    // Note: do not scale DC component
+    psd *= 2;
+    psd.first() /= 2;
+    psd.last() /= 2;
+
+    if (debug)
+    {
+        Pout<< "Average PSD: " << average(psd) << endl;
+    }
+
+    return tpsd;
+}
+
+
+Foam::scalar Foam::noiseModel::RAf(const scalar f) const
+{
+    const scalar c1 = sqr(12194.0);
+    const scalar c2 = sqr(20.6);
+    const scalar c3 = sqr(107.7);
+    const scalar c4 = sqr(739.9);
+
+    const scalar f2 = f*f;
+
+    return
+        c1*sqr(f2)
+       /(
+            (f2 + c2)*sqrt((f2 + c3)*(f2 + c4))*(f2 + c1)
+        );
+}
+
+
+Foam::scalar Foam::noiseModel::gainA(const scalar f) const
+{
+    return 20*log10(RAf(f)) - 20*log10(RAf(1000));
+}
+
+
+Foam::scalar Foam::noiseModel::RBf(const scalar f) const
+{
+    const scalar c1 = sqr(12194.0);
+    const scalar c2 = sqr(20.6);
+    const scalar c3 = sqr(158.5);
+
+    const scalar f2 = f*f;
+
+    return
+        c1*f2*f
+       /(
+            (f2 + c2)*sqrt(f2 + c3)*(f2 + c1)
+        );
+}
+
+
+Foam::scalar Foam::noiseModel::gainB(const scalar f) const
+{
+    return 20*log10(RBf(f)) - 20*log10(RBf(1000));
+}
+
+
+Foam::scalar Foam::noiseModel::RCf(const scalar f) const
+{
+    const scalar c1 = sqr(12194.0);
+    const scalar c2 = sqr(20.6);
+
+    const scalar f2 = f*f;
+
+    return c1*f2/((f2 + c2)*(f2 + c1));
+}
+
+
+Foam::scalar Foam::noiseModel::gainC(const scalar f) const
+{
+    return 20*log10(RCf(f)) - 20*log10(RCf(1000));
+}
+
+
+Foam::scalar Foam::noiseModel::RDf(const scalar f) const
+{
+    const scalar f2 = f*f;
+
+    const scalar hf =
+        (sqr(1037918.48 - f2) + 1080768.16*f2)
+       /(sqr(9837328 - f2) + 11723776*f2);
+
+    return
+        f/6.8966888496476e-5*Foam::sqrt(hf/((f2 + 79919.29)*(f2 + 1345600)));
+}
+
+
+Foam::scalar Foam::noiseModel::gainD(const scalar f) const
+{
+    return 20*log10(RDf(f));
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::noiseModel::noiseModel(const dictionary& dict, const bool readFields)
@@ -164,10 +505,10 @@ Foam::noiseModel::noiseModel(const dictionary& dict, const bool readFields)
     nSamples_(65536),
     fLower_(25),
     fUpper_(10000),
-    customBounds_(false),
     startTime_(0),
     windowModelPtr_(),
     graphFormat_("raw"),
+    SPLweighting_(weightingType::none),
     minPressure_(-0.5*VGREAT),
     maxPressure_(0.5*VGREAT),
     outputPrefix_(),
@@ -175,12 +516,20 @@ Foam::noiseModel::noiseModel(const dictionary& dict, const bool readFields)
     writeSPL_(true),
     writePSD_(true),
     writePSDf_(true),
-    writeOctaves_(true)
+    writeOctaves_(true),
+    planInfo_()
 {
+    planInfo_.active = false;
+
     if (readFields)
     {
         read(dict);
     }
+
+    if (debug)
+    {
+        writeWeightings();
+    }
 }
 
 
@@ -190,15 +539,8 @@ bool Foam::noiseModel::read(const dictionary& dict)
 {
     dict.readIfPresent("rhoRef", rhoRef_);
     dict.readIfPresent("N", nSamples_);
-    customBounds_ = false;
+    dict.readIfPresent("fl", fLower_);
-    if (dict.readIfPresent("fl", fLower_))
+    dict.readIfPresent("fu", fUpper_);
-    {
-        customBounds_ = true;
-    }
-    if (dict.readIfPresent("fu", fUpper_))
-    {
-        customBounds_ = true;
-    }
     dict.readIfPresent("startTime", startTime_);
     dict.readIfPresent("graphFormat", graphFormat_);
     dict.readIfPresent("minPressure", minPressure_);
@@ -210,7 +552,6 @@ bool Foam::noiseModel::read(const dictionary& dict)
         FatalIOErrorInFunction(dict)
             << "fl: lower frequency bound must be greater than zero"
             << exit(FatalIOError);
-
     }
 
     if (fUpper_ < 0)
@@ -218,7 +559,6 @@ bool Foam::noiseModel::read(const dictionary& dict)
         FatalIOErrorInFunction(dict)
             << "fu: upper frequency bound must be greater than zero"
             << exit(FatalIOError);
-
     }
 
     if (fUpper_ < fLower_)
@@ -230,6 +570,14 @@ bool Foam::noiseModel::read(const dictionary& dict)
 
     }
 
+    Info<< "    Frequency bounds:" << nl
+        << "        lower: " << fLower_ << nl
+        << "        upper: " << fUpper_ << endl;
+
+    weightingTypeNames_.readIfPresent("SPLweighting", dict, SPLweighting_);
+
+    Info<< "    Weighting: " << weightingTypeNames_[SPLweighting_] << endl;
+
     Info<< "    Write options:" << endl;
     dictionary optDict(dict.subOrEmptyDict("writeOptions"));
     readWriteOption(optDict, "writePrmsf", writePrmsf_);
@@ -241,10 +589,179 @@ bool Foam::noiseModel::read(const dictionary& dict)
 
     windowModelPtr_ = windowModel::New(dict, nSamples_);
 
+    cleanFFTW();
+
+    const label windowSize = windowModelPtr_->nSamples();
+
+    if (windowSize > 1)
+    {
+        planInfo_.active = true;
+        planInfo_.windowSize = windowSize;
+        planInfo_.in.setSize(windowSize);
+        planInfo_.out.setSize(windowSize);
+
+        // Using real to half-complex fftw 'kind'
+        planInfo_.plan =
+            fftw_plan_r2r_1d
+            (
+                windowSize,
+                planInfo_.in.data(),
+                planInfo_.out.data(),
+                FFTW_R2HC,
+                windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
+            );
+    }
+
     Info<< nl << endl;
 
     return true;
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::noiseModel::PSD
+(
+    const scalarField& PSDf
+) const
+{
+    return 10*log10(PSDf/sqr(2e-5));
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::SPL
+(
+    const scalarField& Prms2,
+    const scalar f
+) const
+{
+    tmp<scalarField> tspl(10*log10(Prms2/sqr(2e-5)));
+    scalarField& spl = tspl.ref();
+
+    switch (SPLweighting_)
+    {
+        case weightingType::none:
+        {
+            break;
+        }
+        case weightingType::dBA:
+        {
+            spl += gainA(f);
+            break;
+        }
+        case weightingType::dBB:
+        {
+            spl += gainB(f);
+            break;
+        }
+        case weightingType::dBC:
+        {
+            spl += gainC(f);
+            break;
+        }
+        case weightingType::dBD:
+        {
+            spl += gainD(f);
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unknown weighting " << weightingTypeNames_[SPLweighting_]
+                << abort(FatalError);
+        }
+    }
+
+    return tspl;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseModel::SPL
+(
+    const scalarField& Prms2,
+    const scalarField& f
+) const
+{
+    tmp<scalarField> tspl(10*log10(Prms2/sqr(2e-5)));
+    scalarField& spl = tspl.ref();
+
+    switch (SPLweighting_)
+    {
+        case weightingType::none:
+        {
+            break;
+        }
+        case weightingType::dBA:
+        {
+            forAll(spl, i)
+            {
+                spl[i] += gainA(f[i]);
+            }
+            break;
+        }
+        case weightingType::dBB:
+        {
+            forAll(spl, i)
+            {
+                spl[i] += gainB(f[i]);
+            }
+            break;
+        }
+        case weightingType::dBC:
+        {
+            forAll(spl, i)
+            {
+                spl[i] += gainC(f[i]);
+            }
+            break;
+        }
+        case weightingType::dBD:
+        {
+            forAll(spl, i)
+            {
+                spl[i] += gainD(f[i]);
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unknown weighting " << weightingTypeNames_[SPLweighting_]
+                << abort(FatalError);
+        }
+    }
+
+    return tspl;
+}
+
+
+void Foam::noiseModel::cleanFFTW()
+{
+    if (planInfo_.active)
+    {
+        planInfo_.active = false;
+        fftw_destroy_plan(planInfo_.plan);
+        fftw_cleanup();
+    }
+}
+
+
+void Foam::noiseModel::writeWeightings() const
+{
+    scalar f0 = 10;
+    scalar f1 = 20000;
+
+    OFstream osA("noiseModel-weight-A");
+    OFstream osB("noiseModel-weight-B");
+    OFstream osC("noiseModel-weight-C");
+    OFstream osD("noiseModel-weight-D");
+
+    for (label f = f0; f <= f1; ++f)
+    {
+        osA << f << " " << gainA(f) << nl;
+        osB << f << " " << gainB(f) << nl;
+        osC << f << " " << gainC(f) << nl;
+        osD << f << " " << gainD(f) << nl;
+    }
+}
+
+
 // ************************************************************************* //

src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2018 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,14 +32,16 @@ Description
     Data is read from a dictionary, e.g.
 
     \verbatim
-    rhoRef          0;
+    rhoRef          1;
     N               4096;
     fl              25;
-    fu              25;
+    fu              10000;
     startTime       0;
 
     outputPrefix    "test1";
 
+    SPLweighting    dBA;
+
     // Optional write options dictionary
     writeOptions
     {
@@ -60,6 +62,7 @@ Description
         fu           | Upper frequency bounds        | no        | 10000
         startTime    | Start time                    | no        | 0
         outputPrefix | Prefix applied to output files| no        | ''
+        SPLweighting | Weighting: dBA, dBB, dBC, DBD | no        | none
         graphFormat  | Graph format                  | no        | raw
         writePrmsf   | Write Prmsf data              | no        | yes
         writeSPL     | Write SPL data                | no        | yes
@@ -80,7 +83,10 @@ SourceFiles
 #include "scalarList.H"
 #include "instantList.H"
 #include "windowModel.H"
+#include "Enum.H"
 #include "runTimeSelectionTables.H"
+#include "graph.H"
+#include <fftw3.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,6 +99,54 @@ namespace Foam
 
 class noiseModel
 {
+public:
+
+    enum class weightingType
+    {
+        none,
+        dBA,
+        dBB,
+        dBC,
+        dBD
+    };
+
+    static const Enum<weightingType> weightingTypeNames_;
+
+    //- FFTW planner information
+    //  Note: storage uses double for use directly with FFTW
+    struct planInfo
+    {
+        bool active;
+        label windowSize;
+        List<double> in;
+        List<double> out;
+        fftw_plan plan;
+    };
+
+    //- Octave band information
+    struct octaveBandInfo
+    {
+        label octave;
+
+        // IDs of bin boundaries in pressure data
+        labelList binIDs;
+
+        // Centre frequencies for each bin
+        scalarField centreFreq;
+    };
+
+
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        noiseModel(const noiseModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const noiseModel&) = delete;
+
+
 protected:
 
     // Protected Data
@@ -112,9 +166,6 @@ protected:
         //- Upper frequency limit, default = 10kHz
         scalar fUpper_;
 
-        //- Flag to indicate that custom frequency bounds are being used
-        bool customBounds_;
-
         //- Start time, default = 0s
         scalar startTime_;
 
@@ -124,6 +175,9 @@ protected:
         //- Graph format
         word graphFormat_;
 
+        //- Weighting
+        weightingType SPLweighting_;
+
 
         // Data validation
 
@@ -155,6 +209,12 @@ protected:
             bool writeOctaves_;
 
 
+        // FFT
+
+            //- Plan information for FFTW
+            mutable planInfo planInfo_;
+
+
     // Protected Member Functions
 
         //- Helper function to read write options and provide info feedback
@@ -184,12 +244,70 @@ protected:
         //- Return the base output directory
         fileName baseFileDir(const label dataseti) const;
 
+        //- Create a field of equally spaced frequencies for the current set of
+        //- data - assumes a constant time step
+        tmp<scalarField> uniformFrequencies(const scalar deltaT) const;
 
-        //- No copy construct
+        //- Return a list of the frequency indices wrt f field that correspond
-        noiseModel(const noiseModel&) = delete;
+        //- to the bands limits for a given octave
+        static void setOctaveBands
+        (
+            const scalarField& f,
+            const scalar fLower,
+            const scalar fUpper,
+            const scalar octave,
+            labelList& fBandIDs,
+            scalarField& fCentre
+        );
 
-        //- No copy assignment
+        //- Generate octave data
-        void operator=(const noiseModel&) = delete;
+        tmp<scalarField> octaves
+        (
+            const scalarField& data,
+            const scalarField& f,
+            const labelUList& freqBandIDs
+        ) const;
+
+        //- Return the fft of the given pressure data
+        tmp<scalarField> Pf(const scalarField& p) const;
+
+        //- Return the multi-window mean fft of the complete pressure data [Pa]
+        tmp<scalarField> meanPf(const scalarField& p) const;
+
+        //- Return the multi-window RMS mean fft of the complete pressure
+        //- data [Pa]
+        tmp<scalarField> RMSmeanPf(const scalarField& p) const;
+
+        //- Return the multi-window Power Spectral Density (PSD) of the complete
+        //- pressure data [Pa^2/Hz]
+        tmp<scalarField> PSDf(const scalarField& p, const scalar deltaT) const;
+
+
+        // Weightings
+
+            //- A weighting function
+            scalar RAf(const scalar f) const;
+
+            //- A weighting as gain in dB
+            scalar gainA(const scalar f) const;
+
+            //- B weighting function
+            scalar RBf(const scalar f) const;
+
+            //- B weighting as gain in dB
+            scalar gainB(const scalar f) const;
+
+            //- C weighting function
+            scalar RCf(const scalar f) const;
+
+            //- C weighting as gain in dB
+            scalar gainC(const scalar f) const;
+
+            //- D weighting function
+            scalar RDf(const scalar f) const;
+
+            //- D weighting as gain in dB
+            scalar gainD(const scalar f) const;
 
 
 public:
@@ -226,6 +344,29 @@ public:
 
         //- Abstract call to calculate
         virtual void calculate() = 0;
+
+        //- PSD [dB/Hz]
+        tmp<Foam::scalarField> PSD(const scalarField& PSDf) const;
+
+        //- SPL [dB]
+        tmp<scalarField> SPL
+        (
+            const scalarField& Prms2,
+            const scalar f
+        ) const;
+
+        //- SPL [dB]
+        tmp<scalarField> SPL
+        (
+            const scalarField& Prms2,
+            const scalarField& f
+        ) const;
+
+        //- Clean up the FFTW
+        void cleanFFTW();
+
+        //- Helper function to check weightings
+        void writeWeightings() const;
 };
 
 

src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2018 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "pointNoise.H"
-#include "noiseFFT.H"
 #include "argList.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -77,12 +76,15 @@ void pointNoise::processData
 {
     Info<< "Reading data file " << data.fName() << endl;
 
-    const word fNameBase = data.fName().nameLessExt();
+    const word fNameBase(data.fName().nameLessExt());
 
     // Time and pressure history data
     scalarField t, p;
     filterTimeData(data.x(), data.y(), t, p);
+
+    // Apply conversions
     p *= rhoRef_;
+    p -= average(p);
 
     Info<< "    read " << t.size() << " values" << nl << endl;
 
@@ -101,79 +103,95 @@ void pointNoise::processData
     windowModelPtr_->validate(t.size());
     const windowModel& win = windowModelPtr_();
     const scalar deltaf = 1.0/(deltaT*win.nSamples());
-    fileName outDir(baseFileDir(dataseti)/fNameBase);
+    const fileName outDir(baseFileDir(dataseti)/fNameBase);
-
-    // Create the fft
-    noiseFFT nfft(deltaT, win.nSamples());
-    nfft.setData(p);
 
 
     // Narrow band data
     // ----------------
 
+    scalarField f(uniformFrequencies(deltaT));
+
     // RMS pressure [Pa]
-    graph Prmsf(nfft.RMSmeanPf(win));
-    if (customBounds_)
-    {
-        Prmsf.setXRange(fLower_, fUpper_);
-    }
     if (writePrmsf_)
     {
-        Info<< "    Creating graph for " << Prmsf.title() << endl;
+        graph g
-        Prmsf.write(outDir, graph::wordify(Prmsf.title()), graphFormat_);
+        (
+            "Prms(f)",
+            "f [Hz]",
+            "Prms(f) [Pa]",
+            f,
+            RMSmeanPf(p)
+        );
+
+        g.setXRange(fLower_, fUpper_);
+
+        Info<< "    Creating graph for " << g.title() << endl;
+        g.write(outDir, graph::wordify(g.title()), graphFormat_);
     }
 
     // PSD [Pa^2/Hz]
-    graph PSDf(nfft.PSDf(win));
+    const scalarField PSDf(this->PSDf(p, deltaT));
-    if (customBounds_)
+
-    {
-        PSDf.setXRange(fLower_, fUpper_);
-    }
     if (writePSDf_)
     {
-        Info<< "    Creating graph for " << PSDf.title() << endl;
+        graph g
-        PSDf.write(outDir, graph::wordify(PSDf.title()), graphFormat_);
+        (
+            "PSD(f)",
+            "f [Hz]",
+            "PSD(f) [PaPa_Hz]",
+            f,
+            PSDf
+        );
+
+        g.setXRange(fLower_, fUpper_);
+
+        Info<< "    Creating graph for " << g.title() << endl;
+        g.write(outDir, graph::wordify(g.title()), graphFormat_);
     }
 
     // PSD [dB/Hz]
-    graph PSDg
-    (
-        "PSD_dB_Hz(f)",
-        "f [Hz]",
-        "PSD(f) [dB_Hz]",
-        Prmsf.x(),
-        noiseFFT::PSD(PSDf.y())
-    );
-
     if (writePSD_)
     {
-        Info<< "    Creating graph for " << PSDg.title() << endl;
+        graph g
-        PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_);
+        (
+            "PSD_dB_Hz(f)",
+            "f [Hz]",
+            "PSD(f) [dB_Hz]",
+            f,
+            PSD(PSDf)
+        );
+
+        g.setXRange(fLower_, fUpper_);
+
+        Info<< "    Creating graph for " << g.title() << endl;
+        g.write(outDir, graph::wordify(g.title()), graphFormat_);
     }
 
     // SPL [dB]
-    graph SPLg
-    (
-        "SPL_dB(f)",
-        "f [Hz]",
-        "SPL(f) [dB]",
-        Prmsf.x(),
-        noiseFFT::SPL(PSDf.y()*deltaf)
-    );
-
     if (writeSPL_)
     {
-        Info<< "    Creating graph for " << SPLg.title() << endl;
+        graph g
-        SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
+        (
+            "SPL_dB(f)",
+            "f [Hz]",
+            "SPL(f) [" + weightingTypeNames_[SPLweighting_] + "]",
+            f,
+            SPL(PSDf*deltaf, f)
+        );
+
+        g.setXRange(fLower_, fUpper_);
+
+        Info<< "    Creating graph for " << g.title() << endl;
+        g.write(outDir, graph::wordify(g.title()), graphFormat_);
     }
 
     if (writeOctaves_)
     {
         labelList octave13BandIDs;
         scalarField octave13FreqCentre;
-        noiseFFT::octaveBandInfo
+        setOctaveBands
         (
-            Prmsf.x(),
+            f,
             fLower_,
             fUpper_,
             3,
@@ -186,18 +204,19 @@ void pointNoise::processData
         // ---------------
 
         // Integrated PSD = P(rms)^2 [Pa^2]
-        graph Prms13f(nfft.octaves(PSDf, octave13BandIDs));
+        scalarField Prms13f(octaves(PSDf, f, octave13BandIDs));
 
-        graph SPL13g
+        graph g
         (
             "SPL13_dB(fm)",
             "fm [Hz]",
-            "SPL(fm) [dB]",
+            "SPL(fm) [" + weightingTypeNames_[SPLweighting_] + "]",
             octave13FreqCentre,
-            noiseFFT::SPL(Prms13f.y())
+            SPL(Prms13f, octave13FreqCentre)
         );
-        Info<< "    Creating graph for " << SPL13g.title() << endl;
+
-        SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_);
+        Info<< "    Creating graph for " << g.title() << endl;
+        g.write(outDir, graph::wordify(g.title()), graphFormat_);
     }
 }
 

src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H

@@ -142,7 +142,6 @@ public:
 
         //- Calculate
         virtual void calculate();
-
 };
 
 

src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C

@@ -28,7 +28,6 @@ License
 #include "surfaceNoise.H"
 #include "surfaceReader.H"
 #include "surfaceWriter.H"
-#include "noiseFFT.H"
 #include "argList.H"
 #include "graph.H"
 #include "addToRunTimeSelectionTable.H"
@@ -192,6 +191,11 @@ void surfaceNoise::readSurfaceData
                 pData[faceI][i] = pSlice[faceI];
             }
         }
+
+        forAll(pData, faceI)
+        {
+            pData[faceI] -= average(pData[faceI]);
+        }
     }
     else
     {
@@ -208,12 +212,21 @@ void surfaceNoise::readSurfaceData
             label timeI = i + startTimeIndex_;
 
             Info<< "    time: " << times_[i] << endl;
-            const scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
+            scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
+
+            // Apply conversions
+            p *= rhoRef_;
+
             forAll(p, faceI)
             {
-                pData[faceI][i] = p[faceI]*rhoRef_;
+                pData[faceI][i] = p[faceI];
             }
         }
+
+        forAll(pData, faceI)
+        {
+            pData[faceI] -= average(pData[faceI]);
+        }
     }
 
     Info<< "Read "
@@ -224,7 +237,7 @@ void surfaceNoise::readSurfaceData
 }
 
 
-Foam::scalar surfaceNoise::writeSurfaceData
+scalar surfaceNoise::writeSurfaceData
 (
     const fileName& outDirBase,
     const word& fName,
@@ -336,7 +349,7 @@ Foam::scalar surfaceNoise::writeSurfaceData
 }
 
 
-Foam::scalar surfaceNoise::surfaceAverage
+scalar surfaceNoise::surfaceAverage
 (
     const scalarField& data,
     const labelList& procFaceOffset
@@ -424,12 +437,6 @@ surfaceNoise::surfaceNoise(const dictionary& dict, const bool readFields)
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-surfaceNoise::~surfaceNoise()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool surfaceNoise::read(const dictionary& dict)
@@ -508,7 +515,7 @@ void surfaceNoise::calculate()
 
         Info<< "Creating noise FFTs" << endl;
 
-        const scalarField freq1(noiseFFT::frequencies(nSamples_, deltaT_));
+        const scalarField freq1(uniformFrequencies(deltaT_));
 
         // Reset desired frequency range if outside actual frequency range
         fLower_ = min(fLower_, max(freq1));
@@ -529,7 +536,7 @@ void surfaceNoise::calculate()
         // Storage for 1/3 octave data
         labelList octave13BandIDs;
         scalarField octave13FreqCentre;
-        noiseFFT::octaveBandInfo
+        setOctaveBands
         (
             freq1,
             fLower_,
@@ -561,32 +568,37 @@ void surfaceNoise::calculate()
         const windowModel& win = windowModelPtr_();
 
         {
-            noiseFFT nfft(deltaT_, win.nSamples());
-
             forAll(pData, faceI)
             {
-                // Note: passes storage from pData to noiseFFT
+                const scalarField& p = pData[faceI];
-                nfft.setData(pData[faceI]);
 
                 // Generate the FFT-based data
-                graph Prmsf(nfft.RMSmeanPf(win));
+                const scalarField Prmsf(RMSmeanPf(p));
-                graph PSDf(nfft.PSDf(win));
+                const scalarField PSDf(this->PSDf(p, deltaT_));
 
                 // Store the frequency results in slot for face of surface
                 forAll(surfPrmsf, i)
                 {
                     label freqI = i*fftWriteInterval_;
-                    surfPrmsf[i][faceI] = Prmsf.y()[freqI];
+                    surfPrmsf[i][faceI] = Prmsf[freqI];
-                    surfPSDf[i][faceI] = PSDf.y()[freqI];
+                    surfPSDf[i][faceI] = PSDf[freqI];
                 }
 
                 // Integrated PSD = P(rms)^2 [Pa^2]
-                graph Prms13f(nfft.octaves(PSDf, octave13BandIDs));
+                const scalarField Prms13f
+                (
+                    octaves
+                    (
+                        PSDf,
+                        freq1,
+                        octave13BandIDs
+                    )
+                );
 
                 // Store the 1/3 octave results in slot for face of surface
                 forAll(surfPrms13f, freqI)
                 {
-                    surfPrms13f[freqI][faceI] = Prms13f.y()[freqI];
+                    surfPrms13f[freqI][faceI] = Prms13f[freqI];
                 }
             }
         }
@@ -663,7 +675,7 @@ void surfaceNoise::calculate()
                         fNameBase,
                         "PSD",
                         freq1[freqI],
-                        noiseFFT::PSD(surfPSDf[i + f0]),
+                        PSD(surfPSDf[i + f0]),
                         procFaceOffset,
                         writePSD_
                     );
@@ -673,7 +685,7 @@ void surfaceNoise::calculate()
                         fNameBase,
                         "SPL",
                         freq1[freqI],
-                        noiseFFT::SPL(surfPSDf[i + f0]*deltaf),
+                        SPL(surfPSDf[i + f0]*deltaf, freq1[freqI]),
                         procFaceOffset,
                         writeSPL_
                     );
@@ -718,7 +730,7 @@ void surfaceNoise::calculate()
                     "f [Hz]",
                     "PSD(f) [dB_Hz]",
                     fOut,
-                    noiseFFT::PSD(PSDfAve)
+                    PSD(PSDfAve)
                 );
                 PSDg.write
                 (
@@ -733,7 +745,7 @@ void surfaceNoise::calculate()
                     "f [Hz]",
                     "SPL(f) [dB]",
                     fOut,
-                    noiseFFT::SPL(PSDfAve*deltaf)
+                    SPL(PSDfAve*deltaf, fOut)
                 );
                 SPLg.write
                 (
@@ -760,7 +772,7 @@ void surfaceNoise::calculate()
                     fNameBase,
                     "SPL13",
                     octave13FreqCentre[i],
-                    noiseFFT::SPL(surfPrms13f[i]),
+                    SPL(surfPrms13f[i], octave13FreqCentre[i]),
                     procFaceOffset,
                     writeOctaves_
                 );
@@ -777,7 +789,7 @@ void surfaceNoise::calculate()
                     "fm [Hz]",
                     "SPL(fm) [dB]",
                     octave13FreqCentre,
-                    noiseFFT::SPL(Prms13fAve)
+                    SPL(Prms13fAve, octave13FreqCentre)
                 );
                 SPL13g.write
                 (

src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H

@@ -211,7 +211,7 @@ public:
     surfaceNoise(const dictionary& dict, const bool readFields = true);
 
     //- Destructor
-    virtual ~surfaceNoise();
+    virtual ~surfaceNoise() = default;
 
 
     // Public Member Functions

src/randomProcesses/windowModels/windowModel/windowModel.H

@@ -30,7 +30,7 @@ Description
     Base class for windowing models
 
 SourceFiles
-    noiseFFT.C
+    windowModel.C
 
 \*---------------------------------------------------------------------------*/