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();
188 scalar deltaT = -1.0;
190 if (times.
size() > 1)
193 deltaT = (times.
last() - times.
first())/scalar(times.
size() - 1);
195 for (label i = 1; i < times.
size(); ++i)
197 scalar dT = times[i] - times[i-1];
199 if (
mag(dT/deltaT - 1) > 1
e-8)
202 <<
"Unable to process data with a variable time step" 210 <<
"Unable to create FFT with a single value" 222 if ((
p[i] < minPressure_) || (
p[i] > maxPressure_))
225 <<
"Pressure data at position " << i
226 <<
" is outside of permitted bounds:" <<
nl 227 <<
" pressure: " <<
p[i] <<
nl 228 <<
" minimum pressure: " << minPressure_ <<
nl 229 <<
" maximum pressure: " << maxPressure_ <<
nl 260 writeHeaderValue(
os,
"Lower frequency", fLower_);
261 writeHeaderValue(
os,
"Upper frequency", fUpper_);
262 writeHeaderValue(
os,
"Window model", windowModelPtr_->type());
263 writeHeaderValue(
os,
"Window number", windowModelPtr_->nWindow());
264 writeHeaderValue(
os,
"Window samples", windowModelPtr_->nSamples());
265 writeHeaderValue(
os,
"Window overlap %", windowModelPtr_->overlapPercent());
266 writeHeaderValue(
os,
"dBRef", dBRef_);
268 for (
const auto& hv : headerValues)
270 writeHeaderValue(
os, hv.first(), hv.second());
273 writeCommented(
os,
x.substr(0,
x.find(
' ')));
274 writeTabbed(
os,
y.substr(0,
y.find(
' ')));
288 if (
f[i] >= fLower_ &&
f[i] <= fUpper_)
302 const auto& window = windowModelPtr_();
303 const label
N = window.nSamples();
308 const scalar deltaf = 1.0/(
N*deltaT);
315 if (
f[i] > fLower_ &&
f[i] < fUpper_)
321 if (
check && nFreq == 0)
324 <<
"No frequencies found in range " 325 << fLower_ <<
" to " << fUpper_
340 if (freqBandIDs.size() < 2)
343 <<
"Octave frequency bands are not defined " 344 <<
"- skipping octaves calculation" 351 auto& octData = toctData.ref();
353 bitSet bandUsed(freqBandIDs.size() - 1);
354 for (label bandI = 0; bandI < freqBandIDs.size() - 1; ++bandI)
356 label fb0 = freqBandIDs[bandI];
357 label fb1 = freqBandIDs[bandI+1];
359 if (fb0 == fb1)
continue;
361 for (label freqI = fb0; freqI < fb1; ++freqI)
364 label f1 =
f[freqI + 1];
365 scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
366 octData[bandI] += dataAve*(f1 - f0);
373 labelList bandUnused = bandUsed.sortedToc();
374 if (bandUnused.size())
377 <<
"Empty bands found: " << bandUnused.size() <<
" of " 387 if (planInfo_.active)
389 if (
p.
size() != planInfo_.windowSize)
392 <<
"Expected pressure data to have " << planInfo_.windowSize
393 <<
" values, but received " <<
p.
size() <<
" values" 397 List<double>& in = planInfo_.in;
398 const List<double>& out = planInfo_.out;
404 ::fftw_execute(planInfo_.plan);
406 const label
n = planInfo_.windowSize;
407 const label nBy2 =
n/2;
409 auto& result = tresult.ref();
414 for (label i = 1; i <= nBy2; ++i)
416 const auto re = out[i];
417 const auto im = out[
n - i];
430 const auto& window = windowModelPtr_();
431 const label
N = window.nSamples();
432 const label nWindow = window.nWindow();
435 auto& meanPf = tmeanPf.ref();
437 for (label windowI = 0; windowI < nWindow; ++windowI)
439 meanPf += Pf(window.apply<scalar>(
p, windowI));
442 meanPf /= scalar(nWindow);
453 const auto& window = windowModelPtr_();
454 const label
N = window.nSamples();
455 const label nWindow = window.nWindow();
458 for (label windowI = 0; windowI < nWindow; ++windowI)
460 RMSMeanPf +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
463 return sqrt(RMSMeanPf/scalar(nWindow))/scalar(
N);
473 const auto& window = windowModelPtr_();
474 const label
N = window.nSamples();
475 const label nWindow = window.nWindow();
478 auto& psd = tpsd.ref();
480 for (label windowI = 0; windowI < nWindow; ++windowI)
482 psd +=
sqr(Pf(window.apply<scalar>(
p, windowI)));
485 scalar fs = 1.0/deltaT;
486 psd /= scalar(nWindow)*fs*
N;
505 const scalar
c1 =
sqr(12194.0);
506 const scalar
c2 =
sqr(20.6);
507 const scalar c3 =
sqr(107.7);
508 const scalar c4 =
sqr(739.9);
510 const scalar f2 =
f*
f;
515 (f2 +
c2)*
sqrt((f2 + c3)*(f2 + c4))*(f2 +
c1)
533 const scalar
c1 =
sqr(12194.0);
534 const scalar
c2 =
sqr(20.6);
535 const scalar c3 =
sqr(158.5);
537 const scalar f2 =
f*
f;
560 const scalar
c1 =
sqr(12194.0);
561 const scalar
c2 =
sqr(20.6);
563 const scalar f2 =
f*
f;
565 return c1*f2/((f2 +
c2)*(f2 +
c1));
582 const scalar f2 =
f*
f;
585 (
sqr(1037918.48 - f2) + 1080768.16*f2)
586 /(
sqr(9837328 - f2) + 11723776*f2);
589 f/6.8966888496476e-5*
Foam::sqrt(hf/((f2 + 79919.29)*(f2 + 1345600)));
614 functionObjects::writeFile(obr,
"noise"),
622 SPLweighting_(weightingType::
none),
624 minPressure_(-0.5*VGREAT),
625 maxPressure_(0.5*VGREAT),
669 <<
"fl: lower frequency bound must be greater than zero" 676 <<
"fu: upper frequency bound must be greater than zero" 680 if (fUpper_ < fLower_)
683 <<
"fu: upper frequency bound must be greater than lower " 684 <<
"frequency bound (fl)" 689 Info<<
" Frequency bounds:" <<
nl 690 <<
" lower: " << fLower_ <<
nl 691 <<
" upper: " << fUpper_ <<
endl;
693 weightingTypeNames_.readIfPresent(
"SPLweighting",
dict, SPLweighting_);
695 Info<<
" Weighting: " << weightingTypeNames_[SPLweighting_] <<
endl;
699 Info<<
" Reference for dB calculation: " << dBRef_ <<
endl;
715 const label windowSize = windowModelPtr_->nSamples();
719 planInfo_.active =
true;
720 planInfo_.windowSize = windowSize;
721 planInfo_.in.setSize(windowSize);
722 planInfo_.out.setSize(windowSize);
730 planInfo_.out.data(),
732 windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
760 switch (SPLweighting_)
766 case weightingType::dBA:
771 case weightingType::dBB:
776 case weightingType::dBC:
781 case weightingType::dBD:
789 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
807 switch (SPLweighting_)
813 case weightingType::dBA:
817 spl[i] += gainA(
f[i]);
821 case weightingType::dBB:
825 spl[i] += gainB(
f[i]);
829 case weightingType::dBC:
833 spl[i] += gainC(
f[i]);
837 case weightingType::dBD:
841 spl[i] += gainD(
f[i]);
848 <<
"Unknown weighting " << weightingTypeNames_[SPLweighting_]
859 if (planInfo_.active)
861 planInfo_.active =
false;
862 fftw_destroy_plan(planInfo_.plan);
873 OFstream osA(
"noiseModel-weight-A");
874 OFstream osB(
"noiseModel-weight-B");
875 OFstream osC(
"noiseModel-weight-C");
876 OFstream osD(
"noiseModel-weight-D");
878 for (label
f = f0;
f <= f1; ++
f)
880 osA <<
f <<
" " << gainA(
f) <<
nl;
881 osB <<
f <<
" " << gainB(
f) <<
nl;
882 osC <<
f <<
" " << gainC(
f) <<
nl;
883 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.
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)
bool none(const UList< bool > &bools)
True if no entries are 'true'.
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)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
T & first()
Access first element of the list, position [0].
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 INVALID.
tmp< scalarField > octaves(const scalarField &data, const scalarField &f, const labelUList &freqBandIDs) const
Generate octave data.
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. ...
T & last()
Access last element of the list, position [size()-1].
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].
messageStream Info
Information stream (stdout output on master, null elsewhere)
void writeWeightings() const
Helper function to check weightings.
tmp< scalarField > safeLog10(const scalarField &fld)
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const wordHashSet &selectedFields, LIFOStack< regIOobject *> &storedObjects)
Read the selected GeometricFields of the templated type.
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)