ENH: Refactored and improved noiseFFT library
This commit is contained in:
andy
parent
02ae553099
commit
7c66e69136
@ -1,25 +1,31 @@
|
||||
Kmesh = Kmesh
|
||||
|
||||
fft = fft
|
||||
|
||||
processes = processes
|
||||
UOprocess = $(processes)/UOprocess
|
||||
|
||||
turbulence = turbulence
|
||||
|
||||
noise = noise
|
||||
|
||||
Kmesh = Kmesh
|
||||
$(Kmesh)/Kmesh.C
|
||||
|
||||
fft = fft
|
||||
$(fft)/fft.C
|
||||
$(fft)/fftRenumber.C
|
||||
$(fft)/calcEk.C
|
||||
$(fft)/kShellIntegration.C
|
||||
|
||||
processes = processes
|
||||
UOprocess = $(processes)/UOprocess
|
||||
$(UOprocess)/UOprocess.C
|
||||
|
||||
turbulence = turbulence
|
||||
$(turbulence)/turbGen.C
|
||||
|
||||
$(noise)/noiseFFT.C
|
||||
noise = noise
|
||||
$(noise)/noiseFFT/noiseFFT.C
|
||||
$(noise)/noiseModels/noiseModel/noiseModel.C
|
||||
$(noise)/noiseModels/noiseModel/noiseModelNew.C
|
||||
$(noise)/noiseModels/pointNoise/pointNoise.C
|
||||
$(noise)/noiseModels/surfaceNoise/surfaceNoise.C
|
||||
|
||||
|
||||
windowModels = windowModels
|
||||
$(windowModels)/windowModel/windowModel.C
|
||||
$(windowModels)/windowModel/windowModelNew.C
|
||||
$(windowModels)/Hanning/Hanning.C
|
||||
|
||||
|
||||
LIB = $(FOAM_LIBBIN)/librandomProcesses
|
||||
|
||||
@ -1,5 +1,9 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/sampling/lnInclude \
|
||||
-I$(LIB_SRC)/surfMesh/lnInclude
|
||||
|
||||
LIB_LIBS = \
|
||||
-lfiniteVolume
|
||||
-lfiniteVolume \
|
||||
-lsampling \
|
||||
-lsurfMesh
|
||||
|
||||
@ -1,454 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "noiseFFT.H"
|
||||
#include "IFstream.H"
|
||||
#include "DynamicList.H"
|
||||
#include "fft.H"
|
||||
#include "SubField.H"
|
||||
#include "mathematicalConstants.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::noiseFFT::p0 = 2e-5;
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::noiseFFT::noiseFFT
|
||||
(
|
||||
const scalar deltat,
|
||||
const scalarField& pressure
|
||||
)
|
||||
:
|
||||
scalarField(pressure),
|
||||
deltat_(deltat)
|
||||
{}
|
||||
|
||||
|
||||
Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
|
||||
:
|
||||
scalarField(),
|
||||
deltat_(0.0)
|
||||
{
|
||||
// Construct pressure data file
|
||||
IFstream pFile(pFileName);
|
||||
|
||||
// Check pFile stream is OK
|
||||
if (!pFile.good())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Cannot read file " << pFileName
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (skip)
|
||||
{
|
||||
scalar dummyt, dummyp;
|
||||
|
||||
for (label i=0; i<skip; i++)
|
||||
{
|
||||
pFile >> dummyt;
|
||||
|
||||
if (!pFile.good() || pFile.eof())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Number of points in file " << pFileName
|
||||
<< " is less than the number to be skipped = " << skip
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
pFile >> dummyp;
|
||||
}
|
||||
}
|
||||
|
||||
scalar t = 0, T = 0;
|
||||
DynamicList<scalar> pData(100000);
|
||||
label i = 0;
|
||||
|
||||
while (!(pFile >> t).eof())
|
||||
{
|
||||
T = t;
|
||||
pFile >> pData(i++);
|
||||
}
|
||||
|
||||
deltat_ = T/pData.size();
|
||||
|
||||
this->transfer(pData);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::graph Foam::noiseFFT::pt() const
|
||||
{
|
||||
scalarField t(size());
|
||||
forAll(t, i)
|
||||
{
|
||||
t[i] = i*deltat_;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"p(t)",
|
||||
"t [s]",
|
||||
"p(t) [Pa]",
|
||||
t,
|
||||
*this
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
|
||||
(
|
||||
const label N,
|
||||
const label ni
|
||||
) const
|
||||
{
|
||||
label windowOffset = N;
|
||||
|
||||
if ((N + ni*windowOffset) > size())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Requested window is outside set of data" << endl
|
||||
<< "number of data = " << size() << endl
|
||||
<< "size of window = " << N << endl
|
||||
<< "window = " << ni
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
tmp<scalarField> tpw(new scalarField(N));
|
||||
scalarField& pw = tpw();
|
||||
|
||||
label offset = ni*windowOffset;
|
||||
|
||||
forAll(pw, i)
|
||||
{
|
||||
pw[i] = operator[](i + offset);
|
||||
}
|
||||
|
||||
return tpw;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
|
||||
{
|
||||
scalarField t(N);
|
||||
forAll(t, i)
|
||||
{
|
||||
t[i] = i*deltat_;
|
||||
}
|
||||
|
||||
scalar T = N*deltat_;
|
||||
|
||||
return 2*(0.5 - 0.5*cos(constant::mathematical::twoPi*t/T));
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
|
||||
(
|
||||
const tmp<scalarField>& tpn
|
||||
) const
|
||||
{
|
||||
tmp<scalarField> tPn2
|
||||
(
|
||||
mag
|
||||
(
|
||||
fft::reverseTransform
|
||||
(
|
||||
ReComplexField(tpn),
|
||||
labelList(1, tpn().size())
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
tpn.clear();
|
||||
|
||||
tmp<scalarField> tPn
|
||||
(
|
||||
new scalarField
|
||||
(
|
||||
scalarField::subField(tPn2(), tPn2().size()/2)
|
||||
)
|
||||
);
|
||||
scalarField& Pn = tPn();
|
||||
|
||||
Pn *= 2.0/sqrt(scalar(tPn2().size()));
|
||||
Pn[0] /= 2.0;
|
||||
|
||||
return tPn;
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::meanPf
|
||||
(
|
||||
const label N,
|
||||
const label nw
|
||||
) const
|
||||
{
|
||||
if (N > size())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Requested window is outside set of data" << nl
|
||||
<< "number of data = " << size() << nl
|
||||
<< "size of window = " << N << nl
|
||||
<< "window = " << nw
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
scalarField MeanPf(N/2, 0.0);
|
||||
|
||||
scalarField Hwf(Hanning(N));
|
||||
|
||||
for (label wi=0; wi<nw; ++wi)
|
||||
{
|
||||
MeanPf += Pf(Hwf*window(N, wi));
|
||||
}
|
||||
|
||||
MeanPf /= nw;
|
||||
|
||||
scalarField f(MeanPf.size());
|
||||
scalar deltaf = 1.0/(N*deltat_);
|
||||
|
||||
forAll(f, i)
|
||||
{
|
||||
f[i] = i*deltaf;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"P(f)",
|
||||
"f [Hz]",
|
||||
"P(f) [Pa]",
|
||||
f,
|
||||
MeanPf
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::RMSmeanPf
|
||||
(
|
||||
const label N,
|
||||
const label nw
|
||||
) const
|
||||
{
|
||||
if (N > size())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Requested window is outside set of data" << endl
|
||||
<< "number of data = " << size() << endl
|
||||
<< "size of window = " << N << endl
|
||||
<< "window = " << nw
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
scalarField RMSMeanPf(N/2, 0.0);
|
||||
|
||||
scalarField Hwf(Hanning(N));
|
||||
|
||||
for (label wi=0; wi<nw; ++wi)
|
||||
{
|
||||
RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
|
||||
}
|
||||
|
||||
RMSMeanPf = sqrt(RMSMeanPf/nw);
|
||||
|
||||
scalarField f(RMSMeanPf.size());
|
||||
scalar deltaf = 1.0/(N*deltat_);
|
||||
|
||||
forAll(f, i)
|
||||
{
|
||||
f[i] = i*deltaf;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"P(f)",
|
||||
"f [Hz]",
|
||||
"P(f) [Pa]",
|
||||
f,
|
||||
RMSMeanPf
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
|
||||
{
|
||||
return graph
|
||||
(
|
||||
"L(f)",
|
||||
"f [Hz]",
|
||||
"L(f) [dB]",
|
||||
gPf.x(),
|
||||
20*log10(gPf.y()/p0)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::Ldelta
|
||||
(
|
||||
const graph& gLf,
|
||||
const scalar f1,
|
||||
const scalar fU
|
||||
) const
|
||||
{
|
||||
const scalarField& f = gLf.x();
|
||||
const scalarField& Lf = gLf.y();
|
||||
|
||||
scalarField ldelta(Lf.size(), 0.0);
|
||||
scalarField fm(ldelta.size());
|
||||
|
||||
scalar fratio = cbrt(2.0);
|
||||
scalar deltaf = 1.0/(2*Lf.size()*deltat_);
|
||||
|
||||
scalar fl = f1/sqrt(fratio);
|
||||
scalar fu = fratio*fl;
|
||||
|
||||
label istart = label(fl/deltaf);
|
||||
label j = 0;
|
||||
|
||||
for (label i = istart; i<Lf.size(); i++)
|
||||
{
|
||||
scalar fmi = sqrt(fu*fl);
|
||||
|
||||
if (fmi > fU + 1) break;
|
||||
|
||||
if (f[i] >= fu)
|
||||
{
|
||||
fm[j] = fmi;
|
||||
ldelta[j] = 10*log10(ldelta[j]);
|
||||
|
||||
j++;
|
||||
|
||||
fl = fu;
|
||||
fu *= fratio;
|
||||
}
|
||||
|
||||
ldelta[j] += pow(10, Lf[i]/10.0);
|
||||
}
|
||||
|
||||
fm.setSize(j);
|
||||
ldelta.setSize(j);
|
||||
|
||||
return graph
|
||||
(
|
||||
"Ldelta",
|
||||
"fm [Hz]",
|
||||
"Ldelta [dB]",
|
||||
fm,
|
||||
ldelta
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::Pdelta
|
||||
(
|
||||
const graph& gPf,
|
||||
const scalar f1,
|
||||
const scalar fU
|
||||
) const
|
||||
{
|
||||
const scalarField& f = gPf.x();
|
||||
const scalarField& Pf = gPf.y();
|
||||
|
||||
scalarField pdelta(Pf.size(), 0.0);
|
||||
scalarField fm(pdelta.size());
|
||||
|
||||
scalar fratio = cbrt(2.0);
|
||||
scalar deltaf = 1.0/(2*Pf.size()*deltat_);
|
||||
|
||||
scalar fl = f1/sqrt(fratio);
|
||||
scalar fu = fratio*fl;
|
||||
|
||||
label istart = label(fl/deltaf + 1.0 - SMALL);
|
||||
label j = 0;
|
||||
|
||||
for (label i = istart; i<Pf.size(); i++)
|
||||
{
|
||||
scalar fmi = sqrt(fu*fl);
|
||||
|
||||
if (fmi > fU + 1) break;
|
||||
|
||||
if (f[i] >= fu)
|
||||
{
|
||||
fm[j] = fmi;
|
||||
pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
|
||||
|
||||
j++;
|
||||
|
||||
fl = fu;
|
||||
fu *= fratio;
|
||||
}
|
||||
|
||||
pdelta[j] += sqr(Pf[i]);
|
||||
}
|
||||
|
||||
fm.setSize(j);
|
||||
pdelta.setSize(j);
|
||||
|
||||
return graph
|
||||
(
|
||||
"Pdelta",
|
||||
"fm [Hz]",
|
||||
"Pdelta [dB]",
|
||||
fm,
|
||||
pdelta
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
|
||||
{
|
||||
const scalarField& Lf = gLf.y();
|
||||
|
||||
scalar lsum = 0.0;
|
||||
|
||||
forAll(Lf, i)
|
||||
{
|
||||
lsum += pow(10, Lf[i]/10.0);
|
||||
}
|
||||
|
||||
lsum = 10*log10(lsum);
|
||||
|
||||
return lsum;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
|
||||
{
|
||||
return p0*pow(10.0, db/20.0);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
|
||||
(
|
||||
const tmp<scalarField>& db
|
||||
) const
|
||||
{
|
||||
return p0*pow(10.0, db/20.0);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
519
src/randomProcesses/noise/noiseFFT/noiseFFT.C
Normal file
519
src/randomProcesses/noise/noiseFFT/noiseFFT.C
Normal file
@ -0,0 +1,519 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
|
||||
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "noiseFFT.H"
|
||||
#include "IFstream.H"
|
||||
#include "DynamicList.H"
|
||||
#include "fft.H"
|
||||
#include "SubField.H"
|
||||
#include "mathematicalConstants.H"
|
||||
|
||||
using namespace Foam::constant;
|
||||
|
||||
// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::noiseFFT::p0 = 2e-5;
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
|
||||
(
|
||||
const label N,
|
||||
const scalar deltaT
|
||||
)
|
||||
{
|
||||
tmp<scalarField> tf(new scalarField(N/2, 0));
|
||||
scalarField& f = tf();
|
||||
|
||||
scalar deltaf = 1.0/(N*deltaT);
|
||||
forAll(f, i)
|
||||
{
|
||||
f[i] = i*deltaf;
|
||||
}
|
||||
|
||||
return tf;
|
||||
}
|
||||
|
||||
|
||||
void Foam::noiseFFT::octaveFrequenciesIDs
|
||||
(
|
||||
const scalarField& f,
|
||||
const scalar fLower,
|
||||
const scalar fUpper,
|
||||
const scalar octave,
|
||||
labelList& freqBandIDs
|
||||
)
|
||||
{
|
||||
// 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^(bandI/octave/2))
|
||||
// Centre frequency given by:
|
||||
// fCentre = f0*(2^(bandI/octave))
|
||||
|
||||
scalar f0 = 1000;
|
||||
scalar minFrequency = max(fLower, min(f));
|
||||
|
||||
// Lower frequency band limit
|
||||
label band0Low = ceil(2*octave*log(minFrequency/f0)/log(2.0));
|
||||
|
||||
// Centre frequency band limit
|
||||
//label band0Centre = ceil(octave*log(fLower/f0)/log(2.0));
|
||||
|
||||
scalar fLowerBand = f0*pow(2, band0Low/octave/2);
|
||||
scalar fRatio = pow(2, 1.0/octave);
|
||||
|
||||
bool complete = false;
|
||||
DynamicList<label> bandIDs(f.size());
|
||||
forAll(f, i)
|
||||
{
|
||||
while (f[i] >= fLowerBand)
|
||||
{
|
||||
bandIDs.append(i);
|
||||
fLowerBand *= fRatio;
|
||||
|
||||
if (fLowerBand > fUpper)
|
||||
{
|
||||
complete = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (complete) break;
|
||||
}
|
||||
|
||||
freqBandIDs.transfer(bandIDs);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::octaveFrequencies
|
||||
(
|
||||
const scalarField& f,
|
||||
const scalar fLower,
|
||||
const scalar fUpper,
|
||||
const scalar octave
|
||||
)
|
||||
{
|
||||
labelList freqBandIDs;
|
||||
octaveFrequenciesIDs(f, fLower, fUpper, octave, freqBandIDs);
|
||||
tmp<scalarField> tf(new scalarField(f, freqBandIDs));
|
||||
return tf;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::noiseFFT::noiseFFT
|
||||
(
|
||||
const scalar deltaT,
|
||||
const scalarField& pressure
|
||||
)
|
||||
:
|
||||
scalarField(pressure),
|
||||
deltaT_(deltaT)
|
||||
{}
|
||||
|
||||
|
||||
Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
|
||||
:
|
||||
scalarField(),
|
||||
deltaT_(0.0)
|
||||
{
|
||||
// Construct pressure data file
|
||||
IFstream pFile(pFileName);
|
||||
|
||||
// Check pFile stream is OK
|
||||
if (!pFile.good())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Cannot read file " << pFileName
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (skip)
|
||||
{
|
||||
scalar dummyt, dummyp;
|
||||
|
||||
for (label i = 0; i < skip; i++)
|
||||
{
|
||||
pFile >> dummyt;
|
||||
|
||||
if (!pFile.good() || pFile.eof())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Number of points in file " << pFileName
|
||||
<< " is less than the number to be skipped = " << skip
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
pFile >> dummyp;
|
||||
}
|
||||
}
|
||||
|
||||
scalar t = 0, T0 = 0, T1 = 0;
|
||||
DynamicList<scalar> pData(100000);
|
||||
label i = 0;
|
||||
|
||||
while (!(pFile >> t).eof())
|
||||
{
|
||||
if (i == 0)
|
||||
{
|
||||
T0 = t;
|
||||
}
|
||||
|
||||
T1 = t;
|
||||
pFile >> pData(i++);
|
||||
}
|
||||
|
||||
// Note: assumes fixed time step
|
||||
deltaT_ = (T1 - T0)/pData.size();
|
||||
|
||||
this->transfer(pData);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::graph Foam::noiseFFT::pt() const
|
||||
{
|
||||
scalarField t(size());
|
||||
forAll(t, i)
|
||||
{
|
||||
t[i] = i*deltaT_;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"p(t)",
|
||||
"t [s]",
|
||||
"p(t) [Pa]",
|
||||
t,
|
||||
*this
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
|
||||
(
|
||||
const tmp<scalarField>& tpn
|
||||
) const
|
||||
{
|
||||
tmp<scalarField> tPn2
|
||||
(
|
||||
mag
|
||||
(
|
||||
fft::reverseTransform
|
||||
(
|
||||
ReComplexField(tpn),
|
||||
labelList(1, tpn().size())
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
tpn.clear();
|
||||
|
||||
tmp<scalarField> tPn
|
||||
(
|
||||
new scalarField
|
||||
(
|
||||
scalarField::subField(tPn2(), tPn2().size()/2)
|
||||
)
|
||||
);
|
||||
scalarField& Pn = tPn();
|
||||
|
||||
Pn *= 2.0/sqrt(scalar(tPn2().size()));
|
||||
Pn[0] /= 2.0;
|
||||
|
||||
return tPn;
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::meanPf(const windowModel& window) const
|
||||
{
|
||||
const label N = window.nSamples();
|
||||
const label nWindow = window.nWindow();
|
||||
|
||||
scalarField meanPf(N/2, 0.0);
|
||||
|
||||
for (label windowI = 0; windowI < nWindow; ++windowI)
|
||||
{
|
||||
meanPf += Pf(window.apply<scalar>(*this, windowI));
|
||||
}
|
||||
|
||||
meanPf /= scalar(nWindow);
|
||||
|
||||
scalar deltaf = 1.0/(N*deltaT_);
|
||||
scalarField f(meanPf.size());
|
||||
forAll(f, i)
|
||||
{
|
||||
f[i] = i*deltaf;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"P(f)",
|
||||
"f [Hz]",
|
||||
"P(f) [Pa]",
|
||||
f,
|
||||
meanPf
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
|
||||
{
|
||||
const label N = window.nSamples();
|
||||
const label nWindow = window.nWindow();
|
||||
|
||||
scalarField RMSMeanPf(N/2, 0.0);
|
||||
|
||||
for (label windowI = 0; windowI < nWindow; ++windowI)
|
||||
{
|
||||
RMSMeanPf += sqr(Pf(window.apply<scalar>(*this, windowI)));
|
||||
}
|
||||
|
||||
RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow));
|
||||
|
||||
scalar deltaf = 1.0/(N*deltaT_);
|
||||
scalarField f(RMSMeanPf.size());
|
||||
forAll(f, i)
|
||||
{
|
||||
f[i] = i*deltaf;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"P(f)",
|
||||
"f [Hz]",
|
||||
"P(f) [Pa]",
|
||||
f,
|
||||
RMSMeanPf
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const
|
||||
{
|
||||
const label N = window.nSamples();
|
||||
const label nWindow = window.nWindow();
|
||||
|
||||
scalarField psd(N/2, 0.0);
|
||||
|
||||
for (label windowI = 0; windowI < nWindow; ++windowI)
|
||||
{
|
||||
psd += 0.5*sqr(Pf(window.apply<scalar>(*this, windowI)));
|
||||
}
|
||||
|
||||
scalar deltaf = 1.0/(N*deltaT_);
|
||||
|
||||
psd /= nWindow*deltaf;
|
||||
|
||||
scalarField f(psd.size());
|
||||
|
||||
forAll(f, i)
|
||||
{
|
||||
f[i] = i*deltaf;
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"PSD(f)",
|
||||
"f [Hz]",
|
||||
"PSD(f) [PaPa_Hz]",
|
||||
f,
|
||||
psd
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::PSD(const graph& gPSD) const
|
||||
{
|
||||
return graph
|
||||
(
|
||||
"PSD(dB)",
|
||||
"f [Hz]",
|
||||
"PSD_dB(f) [dB]",
|
||||
gPSD.x(),
|
||||
10*log10(gPSD.y()/sqr(p0))
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
|
||||
{
|
||||
return graph
|
||||
(
|
||||
"L(f)",
|
||||
"f [Hz]",
|
||||
"L(f) [dB]",
|
||||
gPf.x(),
|
||||
20*log10(gPf.y()/p0)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::Ldelta
|
||||
(
|
||||
const graph& gLf,
|
||||
const labelList& freqBandIDs
|
||||
) const
|
||||
{
|
||||
if (freqBandIDs.size() < 2)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Octave frequency bands are not defined "
|
||||
<< "- skipping Ldelta calculation"
|
||||
<< endl;
|
||||
|
||||
return graph
|
||||
(
|
||||
"Ldelta",
|
||||
"fm [Hz]",
|
||||
"Ldelta [dB]",
|
||||
scalarField(),
|
||||
scalarField()
|
||||
);
|
||||
}
|
||||
|
||||
const scalarField& f = gLf.x();
|
||||
const scalarField& Lf = gLf.y();
|
||||
|
||||
scalarField ldelta(freqBandIDs.size() - 1, 0.0);
|
||||
scalarField fm(freqBandIDs.size() -1, 0.0);
|
||||
|
||||
for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
|
||||
{
|
||||
label f0 = freqBandIDs[bandI];
|
||||
label f1 = freqBandIDs[bandI+1];
|
||||
fm[bandI] = f[f0];
|
||||
|
||||
if (f0 == f1) continue;
|
||||
|
||||
for (label freqI = f0; freqI < f1; freqI++)
|
||||
{
|
||||
ldelta[bandI] += pow(10, Lf[freqI]/10.0);
|
||||
}
|
||||
ldelta[bandI] = 10*log10(ldelta[bandI]);
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"Ldelta",
|
||||
"fm [Hz]",
|
||||
"Ldelta [dB]",
|
||||
fm,
|
||||
ldelta
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::graph Foam::noiseFFT::Pdelta
|
||||
(
|
||||
const graph& gPf,
|
||||
const labelList& freqBandIDs
|
||||
) const
|
||||
{
|
||||
if (freqBandIDs.size() < 2)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Octave frequency bands are not defined "
|
||||
<< "- skipping Pdelta calculation"
|
||||
<< endl;
|
||||
|
||||
return graph
|
||||
(
|
||||
"Pdelta",
|
||||
"fm [Hz]",
|
||||
"Pdelta [dB]",
|
||||
scalarField(),
|
||||
scalarField()
|
||||
);
|
||||
}
|
||||
|
||||
const scalarField& f = gPf.x();
|
||||
const scalarField& Pf = gPf.y();
|
||||
|
||||
scalarField pdelta(freqBandIDs.size() - 1, 0.0);
|
||||
scalarField fm(pdelta.size());
|
||||
|
||||
for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
|
||||
{
|
||||
label f0 = freqBandIDs[bandI];
|
||||
label f1 = freqBandIDs[bandI+1];
|
||||
fm[bandI] = f[f0];
|
||||
|
||||
if (f0 == f1) continue;
|
||||
|
||||
for (label freqI = f0; freqI < f1; freqI++)
|
||||
{
|
||||
pdelta[bandI] += sqr(Pf[freqI]);
|
||||
}
|
||||
pdelta[bandI] = sqrt((2.0/3.0)*pdelta[bandI]);
|
||||
}
|
||||
|
||||
return graph
|
||||
(
|
||||
"Pdelta",
|
||||
"fm [Hz]",
|
||||
"Pdelta [dB]",
|
||||
fm,
|
||||
pdelta
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
|
||||
{
|
||||
const scalarField& Lf = gLf.y();
|
||||
|
||||
scalar lsum = 0.0;
|
||||
|
||||
forAll(Lf, i)
|
||||
{
|
||||
lsum += pow(10, Lf[i]/10.0);
|
||||
}
|
||||
|
||||
lsum = 10*log10(lsum);
|
||||
|
||||
return lsum;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
|
||||
{
|
||||
return p0*pow(10.0, db/20.0);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
|
||||
(
|
||||
const tmp<scalarField>& db
|
||||
) const
|
||||
{
|
||||
return p0*pow(10.0, db/20.0);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -2,8 +2,8 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
|
||||
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -25,11 +25,26 @@ Class
|
||||
Foam::noiseFFT
|
||||
|
||||
Description
|
||||
FFT of the pressure field
|
||||
Performs FFT of pressure field to generate noise data.
|
||||
|
||||
General functionality:
|
||||
- Pf: fft of the pressure data
|
||||
- meanPf: multi-window mean fft
|
||||
- RMSmeanPf: multi-window RMS mean fft
|
||||
- PSDf: multi-window power spectral density (PSD) in frequency domain
|
||||
- PSD: multi-window power spectral density (PSD) in dB
|
||||
- Lf: narrow-band pressure-fluctuation level (PFL) in frequency domain
|
||||
|
||||
Octave-based data:
|
||||
- Ldelta: PFL spectrum
|
||||
- Pdelta: pressure spectrum
|
||||
|
||||
SourceFiles
|
||||
noiseFFT.C
|
||||
|
||||
SeeAlso
|
||||
windowModel.H
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef noiseFFT_H
|
||||
@ -37,6 +52,7 @@ SourceFiles
|
||||
|
||||
#include "scalarField.H"
|
||||
#include "graph.H"
|
||||
#include "windowModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -44,7 +60,7 @@ namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class noiseFFT Declaration
|
||||
Class noiseFFT Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class noiseFFT
|
||||
@ -54,7 +70,7 @@ class noiseFFT
|
||||
// Private data
|
||||
|
||||
//- Time spacing of the raw data
|
||||
scalar deltat_;
|
||||
scalar deltaT_;
|
||||
|
||||
|
||||
public:
|
||||
@ -81,34 +97,72 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Return the FFT frequencies
|
||||
static tmp<scalarField> frequencies
|
||||
(
|
||||
const label N,
|
||||
const scalar deltaT
|
||||
);
|
||||
|
||||
//- Return a list of the frequency indices wrt f field that
|
||||
// correspond to the bands limits for a given octave
|
||||
static void octaveFrequenciesIDs
|
||||
(
|
||||
const scalarField& f,
|
||||
const scalar fLower,
|
||||
const scalar fUpper,
|
||||
const scalar octave,
|
||||
labelList& freqBandIDs
|
||||
);
|
||||
|
||||
//- Return the 1/octave octave frequency bounds
|
||||
static tmp<scalarField> octaveFrequencies
|
||||
(
|
||||
const scalarField& f,
|
||||
const scalar fLower,
|
||||
const scalar fUpper,
|
||||
const scalar octave
|
||||
);
|
||||
|
||||
//- Return the graph of p(t)
|
||||
graph pt() const;
|
||||
|
||||
//- Return the nth window
|
||||
tmp<scalarField> window(const label N, const label n) const;
|
||||
|
||||
//- Return the Hanning window function
|
||||
tmp<scalarField> Hanning(const label N) const;
|
||||
|
||||
//- Return the fft of the given pressure data
|
||||
tmp<scalarField> Pf(const tmp<scalarField>& pn) const;
|
||||
|
||||
//- Return the multi-window mean fft of the complete pressure data
|
||||
graph meanPf(const label N, const label nw) const;
|
||||
//- Return the multi-window mean fft of the complete pressure data [Pa]
|
||||
graph meanPf(const windowModel& window) const;
|
||||
|
||||
//- Return the multi-window RMS mean fft of the complete pressure data
|
||||
graph RMSmeanPf(const label N, const label nw) const;
|
||||
//- Return the multi-window RMS mean fft of the complete pressure
|
||||
// data [Pa]
|
||||
graph RMSmeanPf(const windowModel& window) const;
|
||||
|
||||
//- Return the narrow-band PFL (pressure-fluctuation level) spectrum
|
||||
//- Return the multi-window PSD (power spectral density) of the complete
|
||||
// pressure data [Pa^2/Hz]
|
||||
graph PSDf(const windowModel& window) const;
|
||||
|
||||
//- Return the PSD [dB]
|
||||
graph PSD(const graph& gPSD) const;
|
||||
|
||||
//- Return the narrow-band PFL (pressure-fluctuation level)
|
||||
// spectrum [dB]
|
||||
graph Lf(const graph& gPf) const;
|
||||
|
||||
//- Return the one-third-octave-band PFL spectrum
|
||||
// starting at octave with mean frequency f1
|
||||
graph Ldelta(const graph& gLf, const scalar f1, const scalar fU) const;
|
||||
//- Return the octave-band PFL spectrum starting at octave
|
||||
// frequencies given by the supplied frequency bands [dB]
|
||||
graph Ldelta
|
||||
(
|
||||
const graph& gLf,
|
||||
const labelList& freqBandIDs
|
||||
) const;
|
||||
|
||||
//- Return the one-third-octave-band pressure spectrum
|
||||
// starting at octave with mean frequency f1
|
||||
graph Pdelta(const graph& gLf, const scalar f1, const scalar fU) const;
|
||||
//- Return the octave-band pressure spectrum at octave
|
||||
// frequencies given by the supplied frequency bands [dB]
|
||||
graph Pdelta
|
||||
(
|
||||
const graph& gLf,
|
||||
const labelList& freqBandIDs
|
||||
) const;
|
||||
|
||||
//- Return the total PFL as the sum of Lf over all frequencies
|
||||
scalar Lsum(const graph& gLf) const;
|
||||
Loading…
Reference in New Issue
Block a user