52 const scalar deltaf = 1.0/(
N*deltaT);
91 scalar fTest = 15.625;
93 const scalar fRatio =
pow(2, 1.0/octave);
94 const scalar fRatioL2C =
pow(2, 0.5/octave);
116 if (bandIDs.insert(i))
119 fc.append(fTest*fRatioL2C);
130 fBandIDs = bandIDs.sortedToc();
151 planInfo_.active =
true;
152 planInfo_.windowSize = windowSize;
153 planInfo_.in.setSize(windowSize);
154 planInfo_.out.setSize(windowSize);
162 planInfo_.out.data(),
164 windowSize <= 8192 ? FFTW_MEASURE : FFTW_ESTIMATE
169 planInfo_.active =
false;
178 if (planInfo_.active)
180 planInfo_.active =
false;
181 fftw_destroy_plan(planInfo_.plan);
207 <<
"Cannot read file " << pFileName
213 scalar dummyt, dummyp;
215 for (label i = 0; i < skip; ++i)
219 if (!pFile.good() || pFile.eof())
222 <<
"Number of points in file " << pFileName
223 <<
" is less than the number to be skipped = " << skip
231 scalar t = 0,
T0 = 0, T1 = 0;
232 DynamicList<scalar> pData(1024);
235 while (!(pFile >> t).eof())
247 deltaT_ = (T1 -
T0)/pData.size();
249 this->transfer(pData);
280 if (planInfo_.active)
283 if (pn.
size() != planInfo_.windowSize)
286 <<
"Expected pressure data to have " << planInfo_.windowSize
287 <<
" values, but received " << pn.
size() <<
" values" 291 List<double>& in = planInfo_.in;
292 const List<double>& out = planInfo_.out;
299 ::fftw_execute(planInfo_.plan);
301 const label
n = planInfo_.windowSize;
302 const label nBy2 =
n/2;
304 auto& result = tresult.ref();
309 for (label i = 1; i <= nBy2; ++i)
311 const auto re = out[i];
312 const auto im = out[
n - i];
325 const label
N = window.nSamples();
326 const label nWindow = window.nWindow();
330 for (label windowI = 0; windowI < nWindow; ++windowI)
332 meanPf += Pf(window.apply<scalar>(*
this, windowI));
335 meanPf /= scalar(nWindow);
337 scalar deltaf = 1.0/(
N*deltaT_);
358 const label nWindow = window.
nWindow();
361 for (label windowI = 0; windowI < nWindow; ++windowI)
363 RMSMeanPf +=
sqr(Pf(window.
apply<scalar>(*
this, windowI)));
366 RMSMeanPf =
sqrt(RMSMeanPf/scalar(nWindow))/scalar(
N);
368 scalar deltaf = 1.0/(
N*deltaT_);
389 const label nWindow = window.
nWindow();
393 for (label windowI = 0; windowI < nWindow; ++windowI)
395 psd +=
sqr(Pf(window.
apply<scalar>(*
this, windowI)));
398 scalar deltaf = 1.0/(
N*deltaT_);
399 scalar fs = 1.0/deltaT_;
400 psd /= scalar(nWindow)*fs*
N;
450 if (freqBandIDs.
size() < 2)
453 <<
"Octave frequency bands are not defined " 454 <<
"- skipping octaves calculation" 473 for (label bandI = 0; bandI < freqBandIDs.
size() - 1; ++bandI)
475 label fb0 = freqBandIDs[bandI];
476 label fb1 = freqBandIDs[bandI+1];
479 if (fb0 == fb1)
continue;
481 for (label freqI = fb0; freqI < fb1; ++freqI)
484 label f1 =
f[freqI + 1];
485 scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
486 octData[bandI] += dataAve*(f1 - f0);
503 return p0*
pow(10.0, db/20.0);
512 return p0*
pow(10.0, db/20.0);
Different types of constants.
List< scalar > scalarList
List of scalar.
void size(const label n)
Older name for setAddressableSize.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &f1)
A class for handling file names.
void pop_back(label n=1)
Reduce size by 1 or more elements. Can be called on an empty list.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const scalarField & y() const
graph octaves(const graph &g, const labelUList &freqBandIDs) const
Generate octave data.
const dimensionedScalar re
Classical electron radius: default SI units: [m].
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< scalarField > Pf(const tmp< scalarField > &pn) const
Return the fft of the given pressure data.
label nWindow() const
Return the number of windows.
static tmp< scalarField > frequencies(const label N, const scalar deltaT)
Return the FFT frequencies.
#define forAll(list, i)
Loop across all elements in list.
Class to create, store and output qgraph files.
const dimensionedScalar e
Elementary charge.
~noiseFFT()
Destructor. Cleanup/destroy plan.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static scalar p0
Reference pressure.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
Base class for windowing models.
noiseFFT(const scalar deltaT, const label windowSize=-1)
Construct from pressure field.
label nSamples() const
Return the number of samples in the window.
errorManip< error > abort(error &err)
const uniformDimensionedVectorField & g
graph RMSmeanPf(const windowModel &window) const
Return the multi-window RMS mean fft of the complete pressure.
static tmp< scalarField > SPL(const scalarField &Prms2)
Return the SPL [dB].
Database for solution data, solver performance and other reduced data.
Input from file stream, using an ISstream.
void setData(scalarList &data)
Set the pressure data.
const Vector< label > N(dict.get< Vector< label >>("N"))
static tmp< complexField > realTransform1D(const scalarField &field)
Transform real-value data.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
graph pt() const
Return the graph of pressure as a function of time.
static tmp< scalarField > PSD(const scalarField &PSDf)
Return the PSD [dB/Hz].
graph meanPf(const windowModel &window) const
Return the multi-window mean fft of the complete pressure data [Pa].
tmp< Field< Type > > apply(const Field< Type > &fld, const label windowI) const
Return the windowed data.
#define WarningInFunction
Report a warning using Foam::Warning.
static void octaveBandInfo(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.
void clear() const noexcept
If object pointer points to valid object: delete object and set pointer to nullptr.
A class for managing temporary objects.
dimensionedScalar log10(const dimensionedScalar &ds)
const scalarField & x() const
scalar dbToPa(const scalar db) const
Convert the db into Pa.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const volScalarField & p0
graph PSDf(const windowModel &window) const
Return the multi-window PSD (power spectral density) of the complete.
static constexpr const zero Zero
Global zero (0)