46 {weightingType::none,
"dB"},
47 {weightingType::dBA,
"dBA"},
48 {weightingType::dBB,
"dBB"},
49 {weightingType::dBC,
"dBC"},
50 {weightingType::dBD,
"dBD"},
71 scalar fTest = 15.625;
73 const scalar fRatio =
pow(2, 1.0/octave);
74 const scalar fRatioL2C =
pow(2, 0.5/octave);
80 DynamicList<scalar> fc;
81 DynamicList<scalar> missedBins;
85 while (fTest < fLower)
98 if (stepi) missedBins.append(fTest/fRatio*fRatioL2C);
104 if (bandIDs.insert(i))
107 fc.append(fTest*fRatioL2C);
118 fBandIDs = bandIDs.sortedToc();
120 if (missedBins.size())
122 label nMiss = missedBins.size();
123 label nTotal = nMiss + fc.size() - 1;
125 <<
"Empty bands found: " << nMiss <<
" of " << nTotal
126 <<
" with centre frequencies " <<
flatOutput(missedBins)
135 fCentre.transfer(fc);
145 auto& result = tresult.ref();
190 if (times.
size() > 1)
193 deltaT = (times.
back() - times.
front())/scalar(times.
size() - 1);
195 bool nonUniform =
false;
196 for (label i = 1; i < times.
size() && !nonUniform; ++i)
198 const scalar dT = times[i] - times[i-1];
200 nonUniform = (
mag(dT/deltaT - 1) > 1
e-8);
206 <<
"Detected non-uniform time step:" 207 <<
" resulting FFT may be incorrect" 208 <<
" (or the deltaT is just very small): " << deltaT
215 <<
"Unable to create FFT with 0 or 1 values" 227 if ((
p[i] < minPressure_) || (
p[i] > maxPressure_))
230 <<
"Pressure data at position " << i
231 <<
" is outside of permitted bounds:" <<
nl 232 <<
" pressure: " <<
p[i] <<
nl 233 <<
" minimum pressure: " << minPressure_ <<
nl 234 <<
" maximum pressure: " << maxPressure_ <<
nl 265 writeHeaderValue(
os,
"Lower frequency", fLower_);
266 writeHeaderValue(
os,
"Upper frequency", fUpper_);
267 writeHeaderValue(
os,
"Window model", windowModelPtr_->type());
268 writeHeaderValue(
os,
"Window number", windowModelPtr_->nWindow());
269 writeHeaderValue(
os,
"Window samples", windowModelPtr_->nSamples());
270 writeHeaderValue(
os,
"Window overlap %", windowModelPtr_->overlapPercent());
271 writeHeaderValue(
os,
"dBRef", dBRef_);
273 for (
const auto& hv : headerValues)
275 writeHeaderValue(
os, hv.first(), hv.second());
278 writeCommented(
os,
x.substr(0,
x.find(
' ')));
279 writeTabbed(
os,
y.substr(0,
y.find(
' ')));
293 if (
f[i] >= fLower_ &&
f[i] <= fUpper_)
307 const auto& window = windowModelPtr_();
308 const label
N = window.nSamples();
313 const scalar deltaf = 1.0/(
N*deltaT);
320 if (
f[i] > fLower_ &&
f[i] < fUpper_)
326 if (
check && nFreq == 0)
329 <<
"No frequencies found in range " 330 << fLower_ <<
" to " << fUpper_
345 if (freqBandIDs.size() < 2)
348 <<
"Octave frequency bands are not defined " 349 <<
"- skipping octaves calculation" 356 auto& octData = toctData.ref();
358 bitSet bandUsed(freqBandIDs.size() - 1);
359 for (label bandI = 0; bandI < freqBandIDs.size() - 1; ++bandI)
361 label fb0 = freqBandIDs[bandI];
362 label fb1 = freqBandIDs[bandI+1];
364 if (fb0 == fb1)
continue;
366 for (label freqI = fb0; freqI < fb1; ++freqI)
369 label f1 =
f[freqI + 1];
370 scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
371 octData[bandI] += dataAve*(f1 - f0);
378 labelList bandUnused = bandUsed.sortedToc();
379 if (bandUnused.size())
382 <<
"Empty bands found: " << bandUnused.size() <<
" of " 392 if (planInfo_.active)
394 if (
p.
size() != planInfo_.windowSize)
397 <<
"Expected pressure data to have " << planInfo_.windowSize
398 <<
" values, but received " <<
p.
size() <<
" values" 402 List<double>& in = planInfo_.in;
403 const List<double>& out = planInfo_.out;
409 ::fftw_execute(planInfo_.plan);
411 const label
n = planInfo_.windowSize;
412 const label nBy2 =
n/2;
414 auto& result = tresult.ref();
419 for (label i = 1; i <= nBy2; ++i)
421 const auto re = out[i];
422 const auto im = out[
n - i];
435 const auto& window = windowModelPtr_();
436 const label
N = window.nSamples();
437 const label nWindow = window.nWindow();
440 auto& meanPf = tmeanPf.ref();
442 for (label windowI = 0; windowI < nWindow; ++windowI)
444 meanPf += Pf(window.apply<scalar>(
p, windowI));
447 meanPf /= scalar(nWindow);
458 const auto& window = windowModelPtr_();
459 const label
N = window.nSamples();
460 const label nWindow = window.nWindow();
463 for (label windowI = 0; windowI < nWindow; ++windowI)
465 RMSMeanPf +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
468 return sqrt(RMSMeanPf/scalar(nWindow))/scalar(
N);
478 const auto& window = windowModelPtr_();
479 const label
N = window.nSamples();
480 const label nWindow = window.nWindow();
483 auto& psd = tpsd.ref();
485 for (label windowI = 0; windowI < nWindow; ++windowI)
487 psd +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
490 scalar fs = 1.0/deltaT;
491 psd /= scalar(nWindow)*fs*
N;
510 const scalar
c1 =
sqr(12194.0);
511 const scalar
c2 =
sqr(20.6);
512 const scalar c3 =
sqr(107.7);
513 const scalar c4 =
sqr(739.9);
515 const scalar f2 =
f*
f;
520 (f2 +
c2)*
sqrt((f2 + c3)*(f2 + c4))*(f2 +
c1)
538 const scalar
c1 =
sqr(12194.0);
539 const scalar
c2 =
sqr(20.6);
540 const scalar c3 =
sqr(158.5);
542 const scalar f2 =
f*
f;
565 const scalar
c1 =
sqr(12194.0);
566 const scalar
c2 =
sqr(20.6);
568 const scalar f2 =
f*
f;
570 return c1*f2/((f2 +
c2)*(f2 +
c1));
587 const scalar f2 =
f*
f;
590 (
sqr(1037918.48 - f2) + 1080768.16*f2)
591 /(
sqr(9837328 - f2) + 11723776*f2);
594 f/6.8966888496476e-5*
Foam::sqrt(hf/((f2 + 79919.29)*(f2 + 1345600)));
619 functionObjects::writeFile(obr,
"noise"),
628 SPLweighting_(weightingType::none),
630 minPressure_(-0.5*VGREAT),
631 maxPressure_(0.5*VGREAT),
681 <<
"Lower frequency bound must be greater than zero" 688 <<
"Upper frequency bound must be greater than zero" 692 if (fUpper_ < fLower_)
695 <<
"Upper frequency bound (" << fUpper_
696 <<
") must be greater than lower frequency bound (" 701 Info<<
" Frequency bounds:" <<
nl 702 <<
" lower: " << fLower_ <<
nl 703 <<
" upper: " << fUpper_ <<
nl 715 weightingTypeNames_.readIfPresent(
"SPLweighting",
dict, SPLweighting_);
717 Info<<
" Weighting: " << weightingTypeNames_[SPLweighting_] <<
endl;
721 Info<<
" Reference for dB calculation: " << dBRef_ <<
endl;
737 const label windowSize = windowModelPtr_->nSamples();
741 planInfo_.active =
true;
742 planInfo_.windowSize = windowSize;
743 planInfo_.in.setSize(windowSize);
744 planInfo_.out.setSize(windowSize);
752 planInfo_.out.data(),
754 windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
782 switch (SPLweighting_)
784 case weightingType::none:
788 case weightingType::dBA:
793 case weightingType::dBB:
798 case weightingType::dBC:
803 case weightingType::dBD:
811 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
829 switch (SPLweighting_)
831 case weightingType::none:
835 case weightingType::dBA:
839 spl[i] += gainA(
f[i]);
843 case weightingType::dBB:
847 spl[i] += gainB(
f[i]);
851 case weightingType::dBC:
855 spl[i] += gainC(
f[i]);
859 case weightingType::dBD:
863 spl[i] += gainD(
f[i]);
870 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
881 if (planInfo_.active)
883 planInfo_.active =
false;
884 fftw_destroy_plan(planInfo_.plan);
895 OFstream osA(
"noiseModel-weight-A");
896 OFstream osB(
"noiseModel-weight-B");
897 OFstream osC(
"noiseModel-weight-C");
898 OFstream osD(
"noiseModel-weight-D");
900 for (label
f = f0;
f <= f1; ++
f)
902 osA <<
f <<
" " << gainA(
f) <<
nl;
903 osB <<
f <<
" " << gainB(
f) <<
nl;
904 osC <<
f <<
" " << gainC(
f) <<
nl;
905 osD <<
f <<
" " << gainD(
f) <<
nl;
List< scalar > scalarList
List of scalar.
void cleanFFTW()
Clean up the FFTW.
void size(const label n)
Older name for setAddressableSize.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &f1)
tmp< scalarField > meanPf(const scalarField &p) const
Return the multi-window mean fft of the complete pressure data [Pa].
scalar gainD(const scalar f) const
D weighting as gain in dB.
A class for handling file names.
bool readIfPresentCompat(const word &keyword, std::initializer_list< std::pair< const char *, int >> compat, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val using any compatibility names if needed. FatalIOError if it is found and there are excess tokens.
static void writeHeader(Ostream &os, const word &fieldName)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 2-tuple for storing two objects of dissimilar types. The container is similar in purpose to std::pa...
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionedScalar re
Classical electron radius: default SI units: [m].
constexpr char nl
The newline '\n' character (0x0a)
T & front()
Access first element of the list, position [0].
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
constexpr char tab
The tab '\t' character(0x09)
scalar gainC(const scalar f) const
C weighting as gain in dB.
tmp< scalarField > Pf(const scalarField &p) const
Return the fft of the given pressure data.
Lookup type of boundary radiation properties.
scalar RBf(const scalar f) const
B weighting function.
bool validateBounds(const scalarList &p) const
Return true if all pressure data is within min/max bounds.
UList< label > labelUList
A UList of labels.
#define forAll(list, i)
Loop across all elements in list.
scalar checkUniformTimeStep(const scalarList ×) const
Check and return uniform time step.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
virtual bool read(const dictionary &dict)
Read from dictionary.
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
tmp< scalarField > octaves(const scalarField &data, const scalarField &f, const labelUList &freqBandIDs) const
Generate octave data.
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const NameMatchPredicate &selectedFields, DynamicList< regIOobject *> &storedObjects)
Read the selected GeometricFields of the templated type and store on the objectRegistry.
noiseModel(const noiseModel &)=delete
No copy construct.
void writeFileHeader(Ostream &os, const string &x, const string &y, const UList< Tuple2< string, token >> &headerValues=UList< Tuple2< string, token >>::null()) const
Write output file header.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
scalar RDf(const scalar f) const
D weighting function.
static const Enum< weightingType > weightingTypeNames_
tmp< scalarField > PSDf(const scalarField &p, const scalar deltaT) const
Return the multi-window Power Spectral Density (PSD) of the complete pressure data [Pa^2/Hz]...
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
scalar RAf(const scalar f) const
A weighting function.
void writeFreqDataToFile(Ostream &os, const scalarField &f, const scalarField &fx) const
errorManip< error > abort(error &err)
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
tmp< Foam::scalarField > PSD(const scalarField &PSDf) const
PSD [dB/Hz].
An Ostream is an abstract base class for all output systems (streams, files, token lists...
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
int debug
Static debugging option.
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val. FatalIOError if it is found and the number of tokens i...
planInfo planInfo_
Plan information for FFTW.
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
virtual bool read(const dictionary &dict)
Read.
static void check(const int retVal, const char *what)
dictionary subOrEmptyDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX, const bool mandatory=false) const
Find and return a sub-dictionary as a copy, otherwise return an empty dictionary. ...
const Vector< label > N(dict.get< Vector< label >>("N"))
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static tmp< complexField > realTransform1D(const scalarField &field)
Transform real-value data.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
scalar gainB(const scalar f) const
B weighting as gain in dB.
static void readWriteOption(const dictionary &dict, const word &lookup, bool &option)
fileName baseFileDir() const
Return the base directory for output.
#define WarningInFunction
Report a warning using Foam::Warning.
static autoPtr< windowModel > New(const dictionary &dict, const label nSamples)
Return a reference to the selected window model.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
static void setOctaveBands(const scalarField &f, const scalar fLower, const scalar fUpper, const scalar octave, labelList &fBandIDs, scalarField &fCentre)
Return a list of the frequency indices wrt f field that correspond to the bands limits for a given oc...
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
scalar gainA(const scalar f) const
A weighting as gain in dB.
tmp< scalarField > SPL(const scalarField &Prms2, const scalar f) const
SPL [dB].
T & back()
Access last element of the list, position [size()-1].
messageStream Info
Information stream (stdout output on master, null elsewhere)
void writeWeightings() const
Helper function to check weightings.
tmp< scalarField > safeLog10(const scalarField &fld)
List< label > labelList
A List of labels.
A class for managing temporary objects.
Registry of regIOobjects.
dimensionedScalar log10(const dimensionedScalar &ds)
tmp< scalarField > uniformFrequencies(const scalar deltaT, const bool check) const
Create a field of equally spaced frequencies for the current set of data - assumes a constant time st...
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
tmp< scalarField > RMSmeanPf(const scalarField &p) const
Return the multi-window RMS mean fft of the complete pressure data [Pa].
scalar RCf(const scalar f) const
C weighting function.
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
static constexpr const zero Zero
Global zero (0)