Commit 7972d0b0 authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Noise functionality library and application updates

parent fe13ff56
......@@ -107,7 +107,7 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
word dictName("noiseDict");
fileName dictName(runTime.system()/"noiseDict");
if (args.optionFound("dict"))
{
dictName = args["dict"];
......@@ -118,7 +118,6 @@ int main(int argc, char *argv[])
IOobject
(
dictName,
runTime.system(),
runTime,
IOobject::MUST_READ
)
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object noiseDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
noiseModel surfaceNoise;
surfaceNoiseCoeffs
{
windowModel Hanning;
HanningCoeffs
{
// Window overlap percentage
overlapPercent 50;
symmetric yes;
extended yes;
// Optional number of windows, default = all available
// nWindow 1;
}
/*
windowModel uniform;
uniformCoeffs
{
// Window overlap percentage
overlapPercent 50;
value 1;
// Optional number of windows, default = all available
// nWindow 1;
}
*/
// Input file
inputFile "postProcessing/faceSource1/surface/patch_motorBike_rider-helmet%65/patch_motorBike_rider-helmet%65.case";
// Surface reader
reader ensight;
// Surface writer
writer ensight;
// Collate times for ensight output - ensures geometry is only written once
writeOptions
{
ensight
{
collateTimes 1;
}
}
// Reference density (to convert from kinematic to static pressure)
rhoRef 1.205;
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096; // 8192; // 4096;
// Lower frequency limit, default = 25Hz
//fl 25;
// Upper frequency limit, default = 10kHz
fu 15000;
// Start time, default = 0s
//startTime 0;
// Write interval for FFT data, default = 1
// fftWriteInterval 100;
}
pointNoiseCoeffs
{
csvFileData
{
fileName "pressureData";
nHeaderLine 1;
refColumn 0;
componentColumns (1);
separator " ";
mergeSeparators yes;
}
HanningCoeffs
{
// Window overlap percentage
overlapPercent 50;
symmetric yes;
extended yes;
// Optional number of windows, default = all available
//nWindow 5;
}
// Graph format, default = raw
graphFormat raw;
// Reference density (to convert from kinematic to static pressure)
rhoRef 1.2;
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096;
// Lower frequency limit, default = 25Hz
//fl 25;
// Upper frequency limit, default = 10kHz
//fu 10000;
// Start time, default = 0s
//startTime 0;
// Write interval for FFT data, default = 1
fftWriteInterval 100;
}
// ************************************************************************* //
......@@ -3,7 +3,6 @@ $(Kmesh)/Kmesh.C
fft = fft
$(fft)/fft.C
$(fft)/fftRenumber.C
$(fft)/calcEk.C
$(fft)/kShellIntegration.C
......@@ -26,6 +25,7 @@ windowModels = windowModels
$(windowModels)/windowModel/windowModel.C
$(windowModels)/windowModel/windowModelNew.C
$(windowModels)/Hanning/Hanning.C
$(windowModels)/uniform/uniform.C
LIB = $(FOAM_LIBBIN)/librandomProcesses
EXE_INC = \
-I$(FFTW_ARCH_PATH)/include \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
LIB_LIBS = \
-L$(FFTW_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -lfftw3 \
-lfiniteVolume \
-lsampling \
-lsurfMesh
......@@ -26,9 +26,10 @@ License
#include "noiseFFT.H"
#include "IFstream.H"
#include "DynamicList.H"
#include "fft.H"
#include "SubField.H"
#include "mathematicalConstants.H"
#include "HashSet.H"
#include "fft.H"
using namespace Foam::constant;
......@@ -44,7 +45,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
)
{
tmp<scalarField> tf(new scalarField(N/2, 0));
scalarField& f = tf();
scalarField& f = tf.ref();
scalar deltaf = 1.0/(N*deltaT);
forAll(f, i)
......@@ -56,70 +57,80 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
}
void Foam::noiseFFT::octaveFrequenciesIDs
Foam::tmp<Foam::scalarField> Foam::noiseFFT::PSD(const scalarField& PSDf)
{
return 10*log10(PSDf/sqr(p0));
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::SPL(const scalarField& Prms2)
{
return 10*log10(Prms2/sqr(p0));
}
void Foam::noiseFFT::octaveBandInfo
(
const scalarField& f,
const scalar fLower,
const scalar fUpper,
const scalar octave,
labelList& freqBandIDs
labelList& fBandIDs,
scalarField& fCentre
)
{
// 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))
// fLow = f0*(2^(0.5*bandI/octave))
// Initial (lowest centre frequency)
scalar fTest = 15.625;
scalar f0 = 1000;
scalar minFrequency = max(fLower, min(f));
const scalar fRatio = pow(2, 1.0/octave);
const scalar fRatioL2C = pow(2, 0.5/octave);
// Lower frequency band limit
label band0Low = ceil(2*octave*log(minFrequency/f0)/log(2.0));
// IDs of band IDs
labelHashSet bandIDs(f.size());
// Centre frequency band limit
//label band0Centre = ceil(octave*log(fLower/f0)/log(2.0));
// Centre frequencies
DynamicList<scalar> fc;
scalar fLowerBand = f0*pow(2, band0Low/octave/2);
scalar fRatio = pow(2, 1.0/octave);
// Convert to lower band limit
fTest /= fRatioL2C;
bool complete = false;
DynamicList<label> bandIDs(f.size());
forAll(f, i)
{
while (f[i] >= fLowerBand)
if (f[i] >= fTest)
{
bandIDs.append(i);
fLowerBand *= fRatio;
// Advance band if appropriate
while (f[i] > fTest)
{
fTest *= fRatio;
}
fTest /= fRatio;
if (bandIDs.insert(i))
{
// Also store (next) centre frequency
fc.append(fTest*fRatioL2C);
}
fTest *= fRatio;
if (fLowerBand > fUpper)
if (fTest > fUpper)
{
complete = true;
break;
}
}
if (complete) break;
}
freqBandIDs.transfer(bandIDs);
}
fBandIDs = bandIDs.toc();
// Remove the last centre frequency (beyond upper frequency limit)
fc.remove();
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;
fCentre.transfer(fc);
}
......@@ -133,7 +144,10 @@ Foam::noiseFFT::noiseFFT
:
scalarField(pressure),
deltaT_(deltaT)
{}
{
scalarField& p = *this;
p -= average(p);
}
Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
......@@ -191,6 +205,9 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
deltaT_ = (T1 - T0)/pData.size();
this->transfer(pData);
scalarField& p = *this;
p -= average(p);
}
......@@ -220,6 +237,8 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
const tmp<scalarField>& tpn
) const
{
// Calculate the 2-sided fft
// Note: result not scaled
tmp<scalarField> tPn2
(
mag
......@@ -234,17 +253,15 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
tpn.clear();
// Only storing the positive half
// Note: storing (N/2) values, DC component at position (0)
tmp<scalarField> tPn
(
new scalarField
(
scalarField::subField(tPn2(), tPn2().size()/2)
scalarField::subField(tPn2(), tPn2().size()/2 + 1)
)
);
scalarField& Pn = tPn();
Pn *= 2.0/sqrt(scalar(tPn2().size()));
Pn[0] /= 2.0;
return tPn;
}
......@@ -255,7 +272,7 @@ 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);
scalarField meanPf(N/2 + 1, 0.0);
for (label windowI = 0; windowI < nWindow; ++windowI)
{
......@@ -287,14 +304,13 @@ 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);
scalarField RMSMeanPf(N/2 + 1, 0.0);
for (label windowI = 0; windowI < nWindow; ++windowI)
{
RMSMeanPf += sqr(Pf(window.apply<scalar>(*this, windowI)));
}
RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow));
RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow))/scalar(N);
scalar deltaf = 1.0/(N*deltaT_);
scalarField f(RMSMeanPf.size());
......@@ -305,9 +321,9 @@ Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
return graph
(
"P(f)",
"Prms(f)",
"f [Hz]",
"P(f) [Pa]",
"Prms(f) [Pa]",
f,
RMSMeanPf
);
......@@ -319,19 +335,30 @@ 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);
scalarField psd(N/2 + 1, 0.0);
for (label windowI = 0; windowI < nWindow; ++windowI)
{
psd += 0.5*sqr(Pf(window.apply<scalar>(*this, windowI)));
psd += sqr(Pf(window.apply<scalar>(*this, windowI)));
}
scalar deltaf = 1.0/(N*deltaT_);
scalar fs = 1.0/deltaT_;
psd /= scalar(nWindow)*fs*N;
psd /= nWindow*deltaf;
// Scaling due to use of 1-sided FFT
// Note: do not scale DC component
psd *= 2;
psd.first() /= 2;
psd.last() /= 2;
scalarField f(psd.size());
if (0) // if (debug)
{
Pout<< "Average PSD: " << average(psd) << endl;
}
forAll(f, i)
{
f[i] = i*deltaf;
......@@ -348,159 +375,87 @@ Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const
}
Foam::graph Foam::noiseFFT::PSD(const graph& gPSD) const
Foam::graph Foam::noiseFFT::PSD(const graph& gPSDf) 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)",
"PSD(f)",
"f [Hz]",
"L(f) [dB]",
gPf.x(),
20*log10(gPf.y()/p0)
"PSD_dB(f) [dB_Hz]",
gPSDf.x(),
10*log10(gPSDf.y()/sqr(p0))
);
}
Foam::graph Foam::noiseFFT::Ldelta
Foam::graph Foam::noiseFFT::octaves
(
const graph& gLf,
const labelList& freqBandIDs
const graph& g,
const labelList& freqBandIDs,
bool integrate
) const
{
if (freqBandIDs.size() < 2)
{
WarningInFunction
<< "Octave frequency bands are not defined "
<< "- skipping Ldelta calculation"
<< "- skipping octaves calculation"
<< endl;
return graph
(
"Ldelta",
"fm [Hz]",
"Ldelta [dB]",
"octave",
"x",
"y",
scalarField(),
scalarField()
);
}
const scalarField& f = gLf.x();
const scalarField& Lf = gLf.y();
const scalarField& f = g.x();
const scalarField& data = g.y();
scalarField ldelta(freqBandIDs.size() - 1, 0.0);
scalarField octData(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];
label fb0 = freqBandIDs[bandI];
label fb1 = freqBandIDs[bandI+1];
fm[bandI] = f[fb0];
if (f0 == f1) continue;
if (fb0 == fb1) continue;
for (label freqI = f0; freqI < f1; freqI++)
if (integrate)
{
ldelta[bandI] += pow(10, Lf[freqI]/10.0);
for (label freqI = fb0; freqI < fb1; freqI++)
{
label f0 = f[freqI];
label f1 = f[freqI + 1];
scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
octData[bandI] += dataAve*(f1 - f0);
}
}
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++)
else
{
pdelta[bandI] += sqr(Pf[freqI]);
for (label freqI = fb0; freqI < fb1; freqI++)
{
octData[bandI] += data[freqI];
}
}
pdelta[bandI] = sqrt((2.0/3.0)*pdelta[bandI]);
}
return graph
(
"Pdelta",
"octaves(f)",
"fm [Hz]",
"Pdelta [dB]",
"octave data",
fm,
pdelta
octData
);
}