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