1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
29 #include "DynamicList.H"
31 #include "mathematicalConstants.H"
34 // * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
36 Foam::scalar Foam::noiseFFT::p0 = 2e-5;
38 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 // Construct from pressure field
41 Foam::noiseFFT::noiseFFT
44 const scalarField& pressure
47 scalarField(pressure),
52 // Construct from pressure field file name
53 Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
58 // Construct control dictionary
59 IFstream pFile(pFileName);
61 // Check pFile stream is OK
66 "noiseFFT::noiseFFT(const fileName& pFileName, const label skip)"
67 ) << "Cannot read file " << pFileName
73 scalar dummyt, dummyp;
75 for (label i=0; i<skip; i++)
79 if (!pFile.good() || pFile.eof())
83 "noiseFFT::noiseFFT(const fileName& pFileName, "
85 ) << "Number of points in file " << pFileName
86 << " is less than the number to be skipped = " << skip
95 DynamicList<scalar> pData(100000);
98 while (!(pFile >> t).eof())
104 deltat_ = T/pData.size();
106 scalarField::operator=(pData.shrink());
110 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
112 Foam::graph Foam::noiseFFT::pt() const
114 scalarField t(size());
131 Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
137 label windowOffset = N;
139 if ((N + ni*windowOffset) > size())
141 FatalErrorIn("noiseFFT::window(const label N, const label n) const")
142 << "Requested window is outside set of data" << endl
143 << "number of data = " << size() << endl
144 << "size of window = " << N << endl
149 tmp<scalarField> tpw(new scalarField(N));
150 scalarField& pw = tpw();
152 label offset = ni*windowOffset;
156 pw[i] = operator[](i + offset);
163 Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
171 scalar T = N*deltat_;
173 return 2*(0.5 - 0.5*cos(2*mathematicalConstant::pi*t/T));
177 Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
179 const tmp<scalarField>& tpn
182 tmp<scalarField> tPn2
186 fft::reverseTransform
189 labelList(1, tpn().size())
200 scalarField::subField(tPn2(), tPn2().size()/2)
203 scalarField& Pn = tPn();
205 Pn *= 2.0/sqrt(scalar(tPn2().size()));
212 Foam::graph Foam::noiseFFT::meanPf
220 FatalErrorIn("noiseFFT::meanPf(const label N, const label nw) const")
221 << "Requested window is outside set of data" << endl
222 << "number of data = " << size() << endl
223 << "size of window = " << N << endl
228 scalarField MeanPf(N/2, 0.0);
230 scalarField Hwf = Hanning(N);
232 for (label wi=0; wi<nw; wi++)
234 MeanPf += Pf(Hwf*window(N, wi));
239 scalarField f(MeanPf.size());
240 scalar deltaf = 1.0/(N*deltat_);
258 Foam::graph Foam::noiseFFT::RMSmeanPf
266 FatalErrorIn("noiseFFT::RMSmeanPf(const label N, const label nw) const")
267 << "Requested window is outside set of data" << endl
268 << "number of data = " << size() << endl
269 << "size of window = " << N << endl
274 scalarField RMSMeanPf(N/2, 0.0);
276 scalarField Hwf = Hanning(N);
278 for (label wi=0; wi<nw; wi++)
280 RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
283 RMSMeanPf = sqrt(RMSMeanPf/nw);
285 scalarField f(RMSMeanPf.size());
286 scalar deltaf = 1.0/(N*deltat_);
304 Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
317 Foam::graph Foam::noiseFFT::Ldelta
324 const scalarField& f = gLf.x();
325 const scalarField& Lf = gLf.y();
327 scalarField ldelta(Lf.size(), 0.0);
328 scalarField fm(ldelta.size());
330 scalar fratio = cbrt(2.0);
331 scalar deltaf = 1.0/(2*Lf.size()*deltat_);
333 scalar fl = f1/sqrt(fratio);
334 scalar fu = fratio*fl;
336 label istart = label(fl/deltaf);
339 for (label i = istart; i<Lf.size(); i++)
341 scalar fmi = sqrt(fu*fl);
343 if (fmi > fU + 1) break;
348 ldelta[j] = 10*log10(ldelta[j]);
356 ldelta[j] += pow(10, Lf[i]/10.0);
373 Foam::graph Foam::noiseFFT::Pdelta
380 const scalarField& f = gPf.x();
381 const scalarField& Pf = gPf.y();
383 scalarField pdelta(Pf.size(), 0.0);
384 scalarField fm(pdelta.size());
386 scalar fratio = cbrt(2.0);
387 scalar deltaf = 1.0/(2*Pf.size()*deltat_);
389 scalar fl = f1/sqrt(fratio);
390 scalar fu = fratio*fl;
392 label istart = label(fl/deltaf + 1.0 - SMALL);
395 for (label i = istart; i<Pf.size(); i++)
397 scalar fmi = sqrt(fu*fl);
399 if (fmi > fU + 1) break;
404 pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
412 pdelta[j] += sqr(Pf[i]);
429 Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
431 const scalarField& Lf = gLf.y();
437 lsum += pow(10, Lf[i]/10.0);
440 lsum = 10*log10(lsum);
446 Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
448 return p0*pow(10.0, db/20.0);
452 Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
454 const tmp<scalarField>& db
457 return p0*pow(10.0, db/20.0);
461 // ************************************************************************* //