Go to the documentation of this file.
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();
137 fCentre.transfer(fc);
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);
191 this->transfer(
data);
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;
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"
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];
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);
graph pt() const
Return the graph of pressure as a function of time.
void setData(scalarList &data)
tmp< Field< Type > > apply(const Field< Type > &fld, const label windowI) const
Return the windowed data.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Class to create, store and output qgraph files.
const scalarField & y() const
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.
A class for handling file names.
void clear() const noexcept
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
~noiseFFT()
Destructor. Cleanup/destroy plan.
Input from file stream, using an ISstream.
Different types of constants.
noiseFFT(const scalar deltaT, const label windowSize=-1)
Construct from pressure field.
graph PSDf(const windowModel &window) const
Return the multi-window PSD (power spectral density) of the complete.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const scalarField & x() const
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
scalar dbToPa(const scalar db) const
Convert the db into Pa.
#define forAll(list, i)
Loop across all elements in list.
bool good() const noexcept
True if next operation might succeed.
graph RMSmeanPf(const windowModel &window) const
Return the multi-window RMS mean fft of the complete pressure.
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].
dimensionedScalar log10(const dimensionedScalar &ds)
bool eof() const noexcept
True if end of input seen.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionedScalar re
Classical electron radius: default SI units: [m].
tmp< scalarField > Pf(const tmp< scalarField > &pn) const
Return the fft of the given pressure data.
static tmp< scalarField > frequencies(const label N, const scalar deltaT)
Return the FFT frequencies.
const uniformDimensionedVectorField & g
errorManip< error > abort(error &err)
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static scalar p0
Reference pressure.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label nWindow() const
Return the number of windows.
const dimensionedScalar e
Elementary charge.
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
graph octaves(const graph &g, const labelUList &freqBandIDs) const
Generate octave data.
label nSamples() const
Return the number of samples in the window.
void size(const label n)
Older name for setAddressableSize.
const Vector< label > N(dict.get< Vector< label >>("N"))
const volScalarField & p0
Base class for windowing models.
#define WarningInFunction
Report a warning using Foam::Warning.
Database for solution data, solver performance and other reduced data.
static tmp< complexField > realTransform1D(const scalarField &field)
Transform real-value data.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
static tmp< scalarField > SPL(const scalarField &Prms2)
Return the SPL [dB].