From 7c66e691366a01d86b42a587a4691d5ccd77d477 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 25 Feb 2016 12:55:41 +0000 Subject: [PATCH] ENH: Refactored and improved noiseFFT library --- src/randomProcesses/Make/files | 30 +- src/randomProcesses/Make/options | 8 +- src/randomProcesses/noise/noiseFFT.C | 454 --------------- src/randomProcesses/noise/noiseFFT/noiseFFT.C | 519 ++++++++++++++++++ .../noise/{ => noiseFFT}/noiseFFT.H | 98 +++- 5 files changed, 619 insertions(+), 490 deletions(-) delete mode 100644 src/randomProcesses/noise/noiseFFT.C create mode 100644 src/randomProcesses/noise/noiseFFT/noiseFFT.C rename src/randomProcesses/noise/{ => noiseFFT}/noiseFFT.H (54%) diff --git a/src/randomProcesses/Make/files b/src/randomProcesses/Make/files index 90cb360357..a16bde251c 100644 --- a/src/randomProcesses/Make/files +++ b/src/randomProcesses/Make/files @@ -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 diff --git a/src/randomProcesses/Make/options b/src/randomProcesses/Make/options index 71b7873964..86d0d7baee 100644 --- a/src/randomProcesses/Make/options +++ b/src/randomProcesses/Make/options @@ -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 diff --git a/src/randomProcesses/noise/noiseFFT.C b/src/randomProcesses/noise/noiseFFT.C deleted file mode 100644 index 88cbe3c5ba..0000000000 --- a/src/randomProcesses/noise/noiseFFT.C +++ /dev/null @@ -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 . - -\*---------------------------------------------------------------------------*/ - -#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> 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 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::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 tpw(new scalarField(N)); - scalarField& pw = tpw(); - - label offset = ni*windowOffset; - - forAll(pw, i) - { - pw[i] = operator[](i + offset); - } - - return tpw; -} - - -Foam::tmp 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::noiseFFT::Pf -( - const tmp& tpn -) const -{ - tmp tPn2 - ( - mag - ( - fft::reverseTransform - ( - ReComplexField(tpn), - labelList(1, tpn().size()) - ) - ) - ); - - tpn.clear(); - - tmp 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 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 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 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::noiseFFT::dbToPa -( - const tmp& db -) const -{ - return p0*pow(10.0, db/20.0); -} - - -// ************************************************************************* // diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C new file mode 100644 index 0000000000..0d36c47d16 --- /dev/null +++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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::noiseFFT::frequencies +( + const label N, + const scalar deltaT +) +{ + tmp 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