noiseFFT.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
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
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "noiseFFT.H"
27 #include "IFstream.H"
28 #include "DynamicList.H"
29 #include "fft.H"
30 #include "SubField.H"
31 #include "mathematicalConstants.H"
32 
33 // * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
34 
35 Foam::scalar Foam::noiseFFT::p0 = 2e-5;
36 
37 
38 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
39 
41 (
42  const scalar deltat,
43  const scalarField& pressure
44 )
45 :
46  scalarField(pressure),
47  deltat_(deltat)
48 {}
49 
50 
51 Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
52 :
53  scalarField(),
54  deltat_(0.0)
55 {
56  // Construct pressure data file
57  IFstream pFile(pFileName);
58 
59  // Check pFile stream is OK
60  if (!pFile.good())
61  {
63  << "Cannot read file " << pFileName
64  << exit(FatalError);
65  }
66 
67  if (skip)
68  {
69  scalar dummyt, dummyp;
70 
71  for (label i=0; i<skip; i++)
72  {
73  pFile >> dummyt;
74 
75  if (!pFile.good() || pFile.eof())
76  {
78  << "Number of points in file " << pFileName
79  << " is less than the number to be skipped = " << skip
80  << exit(FatalError);
81  }
82 
83  pFile >> dummyp;
84  }
85  }
86 
87  scalar t = 0, T = 0;
88  DynamicList<scalar> pData(100000);
89  label i = 0;
90 
91  while (!(pFile >> t).eof())
92  {
93  T = t;
94  pFile >> pData(i++);
95  }
96 
97  deltat_ = T/pData.size();
98 
99  this->transfer(pData);
100 }
101 
102 
103 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
104 
106 {
107  scalarField t(size());
108  forAll(t, i)
109  {
110  t[i] = i*deltat_;
111  }
112 
113  return Pair<scalarField>(t, *this);
114 }
115 
116 
118 (
119  const label N,
120  const label ni
121 ) const
122 {
123  label windowOffset = N;
124 
125  if ((N + ni*windowOffset) > size())
126  {
128  << "Requested window is outside set of data" << endl
129  << "number of data = " << size() << endl
130  << "size of window = " << N << endl
131  << "window = " << ni
132  << exit(FatalError);
133  }
134 
135  tmp<scalarField> tpw(new scalarField(N));
136  scalarField& pw = tpw.ref();
137 
138  label offset = ni*windowOffset;
139 
140  forAll(pw, i)
141  {
142  pw[i] = operator[](i + offset);
143  }
144 
145  return tpw;
146 }
147 
148 
150 {
151  scalarField t(N);
152  forAll(t, i)
153  {
154  t[i] = i*deltat_;
155  }
156 
157  scalar T = N*deltat_;
158 
159  return 2*(0.5 - 0.5*cos(constant::mathematical::twoPi*t/T));
160 }
161 
162 
164 (
165  const tmp<scalarField>& tpn
166 ) const
167 {
168  tmp<scalarField> tPn2
169  (
170  mag
171  (
173  (
174  ReComplexField(tpn),
175  labelList(1, tpn().size())
176  )
177  )
178  );
179 
180  tpn.clear();
181 
182  tmp<scalarField> tPn
183  (
184  new scalarField
185  (
186  scalarField::subField(tPn2(), tPn2().size()/2)
187  )
188  );
189  scalarField& Pn = tPn.ref();
190 
191  Pn *= 2.0/sqrt(scalar(tPn2().size()));
192  Pn[0] /= 2.0;
193 
194  return tPn;
195 }
196 
197 
199 (
200  const label N,
201  const label nw
202 ) const
203 {
204  if (N > size())
205  {
207  << "Requested window is outside set of data" << nl
208  << "number of data = " << size() << nl
209  << "size of window = " << N << nl
210  << "window = " << nw
211  << exit(FatalError);
212  }
213 
214  scalarField MeanPf(N/2, 0.0);
215 
216  scalarField Hwf(Hanning(N));
217 
218  for (label wi=0; wi<nw; ++wi)
219  {
220  MeanPf += Pf(Hwf*window(N, wi));
221  }
222 
223  MeanPf /= nw;
224 
225  scalarField f(MeanPf.size());
226  scalar deltaf = 1.0/(N*deltat_);
227 
228  forAll(f, i)
229  {
230  f[i] = i*deltaf;
231  }
232 
233  return Pair<scalarField>(f, MeanPf);
234 }
235 
236 
238 (
239  const label N,
240  const label nw
241 ) const
242 {
243  if (N > size())
244  {
246  << "Requested window is outside set of data" << endl
247  << "number of data = " << size() << endl
248  << "size of window = " << N << endl
249  << "window = " << nw
250  << exit(FatalError);
251  }
252 
253  scalarField RMSMeanPf(N/2, 0.0);
254 
255  scalarField Hwf(Hanning(N));
256 
257  for (label wi=0; wi<nw; ++wi)
258  {
259  RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
260  }
261 
262  RMSMeanPf = sqrt(RMSMeanPf/nw);
263 
264  scalarField f(RMSMeanPf.size());
265  scalar deltaf = 1.0/(N*deltat_);
266 
267  forAll(f, i)
268  {
269  f[i] = i*deltaf;
270  }
271 
272  return Pair<scalarField>(f, RMSMeanPf);
273 }
274 
275 
277 (
278  const Pair<scalarField>& gPf
279 ) const
280 {
281  return Pair<scalarField>(gPf.first(), 20*log10(gPf.second()/p0));
282 }
283 
284 
286 (
287  const Pair<scalarField>& gLf,
288  const scalar f1,
289  const scalar fU
290 ) const
291 {
292  const scalarField& f = gLf.first();
293  const scalarField& Lf = gLf.second();
294 
295  scalarField ldelta(Lf.size(), 0.0);
296  scalarField fm(ldelta.size());
297 
298  scalar fratio = cbrt(2.0);
299  scalar deltaf = 1.0/(2*Lf.size()*deltat_);
300 
301  scalar fl = f1/sqrt(fratio);
302  scalar fu = fratio*fl;
303 
304  label istart = label(fl/deltaf);
305  label j = 0;
306 
307  for (label i = istart; i<Lf.size(); i++)
308  {
309  scalar fmi = sqrt(fu*fl);
310 
311  if (fmi > fU + 1) break;
312 
313  if (f[i] >= fu)
314  {
315  fm[j] = fmi;
316  ldelta[j] = 10*log10(ldelta[j]);
317 
318  j++;
319 
320  fl = fu;
321  fu *= fratio;
322  }
323 
324  ldelta[j] += pow(10, Lf[i]/10.0);
325  }
326 
327  fm.setSize(j);
328  ldelta.setSize(j);
329 
330  return Pair<scalarField>(fm, ldelta);
331 }
332 
333 
335 (
336  const Pair<scalarField>& gPf,
337  const scalar f1,
338  const scalar fU
339 ) const
340 {
341  const scalarField& f = gPf.first();
342  const scalarField& Pf = gPf.second();
343 
344  scalarField pdelta(Pf.size(), 0.0);
345  scalarField fm(pdelta.size());
346 
347  scalar fratio = cbrt(2.0);
348  scalar deltaf = 1.0/(2*Pf.size()*deltat_);
349 
350  scalar fl = f1/sqrt(fratio);
351  scalar fu = fratio*fl;
352 
353  label istart = label(fl/deltaf + 1.0 - small);
354  label j = 0;
355 
356  for (label i = istart; i<Pf.size(); i++)
357  {
358  scalar fmi = sqrt(fu*fl);
359 
360  if (fmi > fU + 1) break;
361 
362  if (f[i] >= fu)
363  {
364  fm[j] = fmi;
365  pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
366 
367  j++;
368 
369  fl = fu;
370  fu *= fratio;
371  }
372 
373  pdelta[j] += sqr(Pf[i]);
374  }
375 
376  fm.setSize(j);
377  pdelta.setSize(j);
378 
379  return Pair<scalarField>(fm, pdelta);
380 }
381 
382 
383 Foam::scalar Foam::noiseFFT::Lsum(const Pair<Foam::scalarField>& gLf) const
384 {
385  const scalarField& Lf = gLf.second();
386 
387  scalar lsum = 0.0;
388 
389  forAll(Lf, i)
390  {
391  lsum += pow(10, Lf[i]/10.0);
392  }
393 
394  lsum = 10*log10(lsum);
395 
396  return lsum;
397 }
398 
399 
400 Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
401 {
402  return p0*pow(10.0, db/20.0);
403 }
404 
405 
407 (
408  const tmp<scalarField>& db
409 ) const
410 {
411  return p0*pow(10.0, db/20.0);
412 }
413 
414 
415 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
tmp< Field< scalar > > T() const
Return the field transpose (only defined for second rank tensors)
Definition: Field.C:539
Input from file stream.
Definition: IFstream.H:85
bool good() const
Return true if next operation might succeed.
Definition: IOstream.H:330
bool eof() const
Return true if end of input seen.
Definition: IOstream.H:336
void transfer(List< scalar > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:65
const Type & second() const
Return second.
Definition: Pair.H:110
const Type & first() const
Return first.
Definition: Pair.H:98
Pre-declare related SubField type.
Definition: SubField.H:63
static tmp< complexField > reverseTransform(const tmp< complexField > &field, const labelList &nn)
Definition: fft.C:203
A class for handling file names.
Definition: fileName.H:82
noiseFFT(const scalar deltat, const scalarField &pressure)
Construct from pressure field.
Definition: noiseFFT.C:41
static scalar p0
Reference pressure.
Definition: noiseFFT.H:61
Pair< scalarField > Ldelta(const Pair< scalarField > &gLf, const scalar f1, const scalar fU) const
Return the one-third-octave-band PFL spectrum.
Definition: noiseFFT.C:286
Pair< scalarField > pt() const
Return p(t)
Definition: noiseFFT.C:105
Pair< scalarField > RMSmeanPf(const label N, const label nw) const
Return the multi-window RMS mean fft of the complete pressure data.
Definition: noiseFFT.C:238
scalar Lsum(const Pair< scalarField > &gLf) const
Return the total PFL as the sum of Lf over all frequencies.
Definition: noiseFFT.C:383
Pair< scalarField > Pdelta(const Pair< scalarField > &gLf, const scalar f1, const scalar fU) const
Return the one-third-octave-band pressure spectrum.
Definition: noiseFFT.C:335
Pair< scalarField > meanPf(const label N, const label nw) const
Return the multi-window mean fft of the complete pressure data.
Definition: noiseFFT.C:199
tmp< scalarField > Hanning(const label N) const
Return the Hanning window function.
Definition: noiseFFT.C:149
scalar dbToPa(const scalar db) const
Convert the db into Pa.
Definition: noiseFFT.C:400
Pair< scalarField > Lf(const Pair< scalarField > &gPf) const
Return the narrow-band PFL (pressure-fluctuation level) spectrum.
Definition: noiseFFT.C:277
tmp< scalarField > window(const label N, const label n) const
Return the nth window.
Definition: noiseFFT.C:118
tmp< scalarField > Pf(const tmp< scalarField > &pn) const
Return the fft of the given pressure data.
Definition: noiseFFT.C:164
A class for managing temporary objects.
Definition: tmp.H:55
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const dimensionedScalar e
Elementary charge.
const scalar twoPi(2 *pi)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
List< label > labelList
A List of labels.
Definition: labelList.H:56
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
dimensionedScalar log10(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
complexField ReComplexField(const UList< scalar > &sf)
Definition: complexFields.C:56
dimensioned< scalar > mag(const dimensioned< Type > &)
error FatalError
dimensionedScalar cbrt(const dimensionedScalar &ds)
void offset(label &lst, const label o)
static const char nl
Definition: Ostream.H:266
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensionedScalar cos(const dimensionedScalar &ds)
labelList f(nPoints)