From 67204543d0ddd9dc0efba0d742e3f8594812e405 Mon Sep 17 00:00:00 2001 From: Andrew Heather <> Date: Tue, 8 Dec 2020 17:35:36 +0000 Subject: [PATCH] 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. --- .../postProcessing/noise/Make/options | 3 +- .../utilities/postProcessing/noise/noise.C | 1 - src/randomProcesses/noise/noiseFFT/noiseFFT.C | 2 +- .../noise/noiseModels/noiseModel/noiseModel.C | 557 +++++++++++++++++- .../noise/noiseModels/noiseModel/noiseModel.H | 161 ++++- .../noise/noiseModels/pointNoise/pointNoise.C | 123 ++-- .../noise/noiseModels/pointNoise/pointNoise.H | 1 - .../noiseModels/surfaceNoise/surfaceNoise.C | 70 ++- .../noiseModels/surfaceNoise/surfaceNoise.H | 2 +- .../windowModels/windowModel/windowModel.H | 2 +- 10 files changed, 805 insertions(+), 117 deletions(-) diff --git a/applications/utilities/postProcessing/noise/Make/options b/applications/utilities/postProcessing/noise/Make/options index dfbcbf2403..92f606cf1a 100644 --- a/applications/utilities/postProcessing/noise/Make/options +++ b/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 \ diff --git a/applications/utilities/postProcessing/noise/noise.C b/applications/utilities/postProcessing/noise/noise.C index ee4681b5f5..ebab1811bc 100644 --- a/applications/utilities/postProcessing/noise/noise.C +++ b/applications/utilities/postProcessing/noise/noise.C @@ -87,7 +87,6 @@ Usage See also - noiseFFT.H noiseModel.H windowModel.H diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C index a2cfbd915a..946a54e79e 100644 --- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C +++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C @@ -127,7 +127,7 @@ void Foam::noiseFFT::octaveBandInfo } } - fBandIDs = bandIDs.toc(); + fBandIDs = bandIDs.sortedToc(); if (fc.size()) { diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C index 7c52fcdb86..82cd1e9c4f 100644 --- a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C +++ b/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::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 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::noiseModel::uniformFrequencies +( + const scalar deltaT +) const +{ + const auto& window = windowModelPtr_(); + const label N = window.nSamples(); + + auto tf = tmp::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::noiseModel::octaves +( + const scalarField& data, + const scalarField& f, + const labelUList& freqBandIDs +) const +{ + auto toctData = tmp::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::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& in = planInfo_.in; + const List& 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::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::noiseModel::meanPf(const scalarField& p) const +{ + const auto& window = windowModelPtr_(); + const label N = window.nSamples(); + const label nWindow = window.nWindow(); + + auto tmeanPf = tmp::New(N/2 + 1, Zero); + auto& meanPf = tmeanPf.ref(); + + for (label windowI = 0; windowI < nWindow; ++windowI) + { + meanPf += Pf(window.apply(p, windowI)); + } + + meanPf /= scalar(nWindow); + + return tmeanPf; +} + + +Foam::tmp 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(p, windowI))); + } + + return sqrt(RMSMeanPf/scalar(nWindow))/scalar(N); +} + + +Foam::tmp 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::New(N/2 + 1, Zero); + auto& psd = tpsd.ref(); + + for (label windowI = 0; windowI < nWindow; ++windowI) + { + psd += sqr(Pf(window.apply(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; - if (dict.readIfPresent("fl", fLower_)) - { - customBounds_ = true; - } - if (dict.readIfPresent("fu", fUpper_)) - { - customBounds_ = true; - } + dict.readIfPresent("fl", fLower_); + dict.readIfPresent("fu", fUpper_); 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::noiseModel::PSD +( + const scalarField& PSDf +) const +{ + return 10*log10(PSDf/sqr(2e-5)); +} + + +Foam::tmp Foam::noiseModel::SPL +( + const scalarField& Prms2, + const scalar f +) const +{ + tmp 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::noiseModel::SPL +( + const scalarField& Prms2, + const scalarField& f +) const +{ + tmp 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; + } +} + + // ************************************************************************* // diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H index e04d137d71..380ec389d2 100644 --- a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H +++ b/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 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -93,6 +99,54 @@ namespace Foam class noiseModel { +public: + + enum class weightingType + { + none, + dBA, + dBB, + dBC, + dBD + }; + + static const Enum weightingTypeNames_; + + //- FFTW planner information + // Note: storage uses double for use directly with FFTW + struct planInfo + { + bool active; + label windowSize; + List in; + List 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 uniformFrequencies(const scalar deltaT) const; - //- No copy construct - noiseModel(const noiseModel&) = delete; + //- Return a list of the frequency indices wrt f field that correspond + //- 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 - void operator=(const noiseModel&) = delete; + //- Generate octave data + tmp octaves + ( + const scalarField& data, + const scalarField& f, + const labelUList& freqBandIDs + ) const; + + //- Return the fft of the given pressure data + tmp Pf(const scalarField& p) const; + + //- Return the multi-window mean fft of the complete pressure data [Pa] + tmp meanPf(const scalarField& p) const; + + //- Return the multi-window RMS mean fft of the complete pressure + //- data [Pa] + tmp RMSmeanPf(const scalarField& p) const; + + //- Return the multi-window Power Spectral Density (PSD) of the complete + //- pressure data [Pa^2/Hz] + tmp 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 PSD(const scalarField& PSDf) const; + + //- SPL [dB] + tmp SPL + ( + const scalarField& Prms2, + const scalar f + ) const; + + //- SPL [dB] + tmp SPL + ( + const scalarField& Prms2, + const scalarField& f + ) const; + + //- Clean up the FFTW + void cleanFFTW(); + + //- Helper function to check weightings + void writeWeightings() const; }; diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C index 80a7a02ff2..8606462f1e 100644 --- a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C +++ b/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); - - // Create the fft - noiseFFT nfft(deltaT, win.nSamples()); - nfft.setData(p); + const fileName outDir(baseFileDir(dataseti)/fNameBase); // 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; - Prmsf.write(outDir, graph::wordify(Prmsf.title()), graphFormat_); + graph g + ( + "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)); - if (customBounds_) - { - PSDf.setXRange(fLower_, fUpper_); - } + const scalarField PSDf(this->PSDf(p, deltaT)); + if (writePSDf_) { - Info<< " Creating graph for " << PSDf.title() << endl; - PSDf.write(outDir, graph::wordify(PSDf.title()), graphFormat_); + graph g + ( + "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; - PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_); + graph g + ( + "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; - SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_); + graph g + ( + "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_); } } diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H index 0296ac2b71..99d2cf8629 100644 --- a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H +++ b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H @@ -142,7 +142,6 @@ public: //- Calculate virtual void calculate(); - }; diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C index a3d4bd79bb..ea227417c4 100644 --- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C +++ b/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 - nfft.setData(pData[faceI]); + const scalarField& p = pData[faceI]; // Generate the FFT-based data - graph Prmsf(nfft.RMSmeanPf(win)); - graph PSDf(nfft.PSDf(win)); + const scalarField Prmsf(RMSmeanPf(p)); + 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]; - surfPSDf[i][faceI] = PSDf.y()[freqI]; + surfPrmsf[i][faceI] = Prmsf[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 ( diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H index 7854839b02..17a589cc29 100644 --- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H +++ b/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 diff --git a/src/randomProcesses/windowModels/windowModel/windowModel.H b/src/randomProcesses/windowModels/windowModel/windowModel.H index c4cb5aec46..22bb440338 100644 --- a/src/randomProcesses/windowModels/windowModel/windowModel.H +++ b/src/randomProcesses/windowModels/windowModel/windowModel.H @@ -30,7 +30,7 @@ Description Base class for windowing models SourceFiles - noiseFFT.C + windowModel.C \*---------------------------------------------------------------------------*/