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-2018 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 graph
114  (
115  "p(t)",
116  "t [s]",
117  "p(t) [Pa]",
118  t,
119  *this
120  );
121 }
122 
123 
125 (
126  const label N,
127  const label ni
128 ) const
129 {
130  label windowOffset = N;
131 
132  if ((N + ni*windowOffset) > size())
133  {
135  << "Requested window is outside set of data" << endl
136  << "number of data = " << size() << endl
137  << "size of window = " << N << endl
138  << "window = " << ni
139  << exit(FatalError);
140  }
141 
142  tmp<scalarField> tpw(new scalarField(N));
143  scalarField& pw = tpw.ref();
144 
145  label offset = ni*windowOffset;
146 
147  forAll(pw, i)
148  {
149  pw[i] = operator[](i + offset);
150  }
151 
152  return tpw;
153 }
154 
155 
157 {
158  scalarField t(N);
159  forAll(t, i)
160  {
161  t[i] = i*deltat_;
162  }
163 
164  scalar T = N*deltat_;
165 
166  return 2*(0.5 - 0.5*cos(constant::mathematical::twoPi*t/T));
167 }
168 
169 
171 (
172  const tmp<scalarField>& tpn
173 ) const
174 {
175  tmp<scalarField> tPn2
176  (
177  mag
178  (
180  (
181  ReComplexField(tpn),
182  labelList(1, tpn().size())
183  )
184  )
185  );
186 
187  tpn.clear();
188 
189  tmp<scalarField> tPn
190  (
191  new scalarField
192  (
193  scalarField::subField(tPn2(), tPn2().size()/2)
194  )
195  );
196  scalarField& Pn = tPn.ref();
197 
198  Pn *= 2.0/sqrt(scalar(tPn2().size()));
199  Pn[0] /= 2.0;
200 
201  return tPn;
202 }
203 
204 
206 (
207  const label N,
208  const label nw
209 ) const
210 {
211  if (N > size())
212  {
214  << "Requested window is outside set of data" << nl
215  << "number of data = " << size() << nl
216  << "size of window = " << N << nl
217  << "window = " << nw
218  << exit(FatalError);
219  }
220 
221  scalarField MeanPf(N/2, 0.0);
222 
223  scalarField Hwf(Hanning(N));
224 
225  for (label wi=0; wi<nw; ++wi)
226  {
227  MeanPf += Pf(Hwf*window(N, wi));
228  }
229 
230  MeanPf /= nw;
231 
232  scalarField f(MeanPf.size());
233  scalar deltaf = 1.0/(N*deltat_);
234 
235  forAll(f, i)
236  {
237  f[i] = i*deltaf;
238  }
239 
240  return graph
241  (
242  "P(f)",
243  "f [Hz]",
244  "P(f) [Pa]",
245  f,
246  MeanPf
247  );
248 }
249 
250 
252 (
253  const label N,
254  const label nw
255 ) const
256 {
257  if (N > size())
258  {
260  << "Requested window is outside set of data" << endl
261  << "number of data = " << size() << endl
262  << "size of window = " << N << endl
263  << "window = " << nw
264  << exit(FatalError);
265  }
266 
267  scalarField RMSMeanPf(N/2, 0.0);
268 
269  scalarField Hwf(Hanning(N));
270 
271  for (label wi=0; wi<nw; ++wi)
272  {
273  RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
274  }
275 
276  RMSMeanPf = sqrt(RMSMeanPf/nw);
277 
278  scalarField f(RMSMeanPf.size());
279  scalar deltaf = 1.0/(N*deltat_);
280 
281  forAll(f, i)
282  {
283  f[i] = i*deltaf;
284  }
285 
286  return graph
287  (
288  "P(f)",
289  "f [Hz]",
290  "P(f) [Pa]",
291  f,
292  RMSMeanPf
293  );
294 }
295 
296 
298 {
299  return graph
300  (
301  "L(f)",
302  "f [Hz]",
303  "L(f) [dB]",
304  gPf.x(),
305  20*log10(gPf.y()/p0)
306  );
307 }
308 
309 
311 (
312  const graph& gLf,
313  const scalar f1,
314  const scalar fU
315 ) const
316 {
317  const scalarField& f = gLf.x();
318  const scalarField& Lf = gLf.y();
319 
320  scalarField ldelta(Lf.size(), 0.0);
321  scalarField fm(ldelta.size());
322 
323  scalar fratio = cbrt(2.0);
324  scalar deltaf = 1.0/(2*Lf.size()*deltat_);
325 
326  scalar fl = f1/sqrt(fratio);
327  scalar fu = fratio*fl;
328 
329  label istart = label(fl/deltaf);
330  label j = 0;
331 
332  for (label i = istart; i<Lf.size(); i++)
333  {
334  scalar fmi = sqrt(fu*fl);
335 
336  if (fmi > fU + 1) break;
337 
338  if (f[i] >= fu)
339  {
340  fm[j] = fmi;
341  ldelta[j] = 10*log10(ldelta[j]);
342 
343  j++;
344 
345  fl = fu;
346  fu *= fratio;
347  }
348 
349  ldelta[j] += pow(10, Lf[i]/10.0);
350  }
351 
352  fm.setSize(j);
353  ldelta.setSize(j);
354 
355  return graph
356  (
357  "Ldelta",
358  "fm [Hz]",
359  "Ldelta [dB]",
360  fm,
361  ldelta
362  );
363 }
364 
365 
367 (
368  const graph& gPf,
369  const scalar f1,
370  const scalar fU
371 ) const
372 {
373  const scalarField& f = gPf.x();
374  const scalarField& Pf = gPf.y();
375 
376  scalarField pdelta(Pf.size(), 0.0);
377  scalarField fm(pdelta.size());
378 
379  scalar fratio = cbrt(2.0);
380  scalar deltaf = 1.0/(2*Pf.size()*deltat_);
381 
382  scalar fl = f1/sqrt(fratio);
383  scalar fu = fratio*fl;
384 
385  label istart = label(fl/deltaf + 1.0 - small);
386  label j = 0;
387 
388  for (label i = istart; i<Pf.size(); i++)
389  {
390  scalar fmi = sqrt(fu*fl);
391 
392  if (fmi > fU + 1) break;
393 
394  if (f[i] >= fu)
395  {
396  fm[j] = fmi;
397  pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
398 
399  j++;
400 
401  fl = fu;
402  fu *= fratio;
403  }
404 
405  pdelta[j] += sqr(Pf[i]);
406  }
407 
408  fm.setSize(j);
409  pdelta.setSize(j);
410 
411  return graph
412  (
413  "Pdelta",
414  "fm [Hz]",
415  "Pdelta [dB]",
416  fm,
417  pdelta
418  );
419 }
420 
421 
422 Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
423 {
424  const scalarField& Lf = gLf.y();
425 
426  scalar lsum = 0.0;
427 
428  forAll(Lf, i)
429  {
430  lsum += pow(10, Lf[i]/10.0);
431  }
432 
433  lsum = 10*log10(lsum);
434 
435  return lsum;
436 }
437 
438 
439 Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
440 {
441  return p0*pow(10.0, db/20.0);
442 }
443 
444 
446 (
447  const tmp<scalarField>& db
448 ) const
449 {
450  return p0*pow(10.0, db/20.0);
451 }
452 
453 
454 // ************************************************************************* //
noiseFFT(const scalar deltat, const scalarField &pressure)
Construct from pressure field.
Definition: noiseFFT.C:41
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
A class for handling file names.
Definition: fileName.H:79
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
scalar & operator[](const label)
Return element of UList.
Definition: UListI.H:167
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
graph Pdelta(const graph &gLf, const scalar f1, const scalar fU) const
Return the one-third-octave-band pressure spectrum.
Definition: noiseFFT.C:367
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const scalarField & y() const
Definition: graph.C:137
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
tmp< scalarField > Pf(const tmp< scalarField > &pn) const
Return the fft of the given pressure data.
Definition: noiseFFT.C:171
static tmp< complexField > reverseTransform(const tmp< complexField > &field, const labelList &nn)
Definition: fft.C:203
Pre-declare related SubField type.
Definition: Field.H:60
graph meanPf(const label N, const label nw) const
Return the multi-window mean fft of the complete pressure data.
Definition: noiseFFT.C:206
bool good() const
Return true if next operation might succeed.
Definition: IOstream.H:330
Class to create, store and output qgraph files.
Definition: graph.H:58
tmp< scalarField > Hanning(const label N) const
Return the Hanning window function.
Definition: noiseFFT.C:156
graph RMSmeanPf(const label N, const label nw) const
Return the multi-window RMS mean fft of the complete pressure data.
Definition: noiseFFT.C:252
dimensionedScalar cos(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static scalar p0
Reference pressure.
Definition: noiseFFT.H:62
dimensionedScalar cbrt(const dimensionedScalar &ds)
graph Ldelta(const graph &gLf, const scalar f1, const scalar fU) const
Return the one-third-octave-band PFL spectrum.
Definition: noiseFFT.C:311
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool eof() const
Return true if end of input seen.
Definition: IOstream.H:336
const scalar twoPi(2 *pi)
volScalarField scalarField(fieldObject, mesh)
static const char nl
Definition: Ostream.H:260
Input from file stream.
Definition: IFstream.H:81
labelList f(nPoints)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
complexField ReComplexField(const UList< scalar > &sf)
Definition: complexFields.C:56
void setSize(const label)
Reset size of List.
Definition: List.C:281
graph pt() const
Return the graph of p(t)
Definition: noiseFFT.C:105
scalar Lsum(const graph &gLf) const
Return the total PFL as the sum of Lf over all frequencies.
Definition: noiseFFT.C:422
tmp< scalarField > window(const label N, const label n) const
Return the nth window.
Definition: noiseFFT.C:125
tmp< Field< scalar > > T() const
Return the field transpose (only defined for second rank tensors)
Definition: Field.C:501
label N
Definition: createFields.H:9
dimensioned< scalar > mag(const dimensioned< Type > &)
A class for managing temporary objects.
Definition: PtrList.H:53
label size() const
Return the number of elements in the UList.
Definition: ListI.H:171
dimensionedScalar log10(const dimensionedScalar &ds)
graph Lf(const graph &gPf) const
Return the narrow-band PFL (pressure-fluctuation level) spectrum.
Definition: noiseFFT.C:297
void transfer(List< scalar > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
const scalarField & x() const
Definition: graph.H:165
scalar dbToPa(const scalar db) const
Convert the db into Pa.
Definition: noiseFFT.C:439
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
label nw
Definition: createFields.H:12