FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/IFstream.H>
28 #include <OpenFOAM/DynamicList.H>
29 #include <randomProcesses/fft.H>
31 #include <OpenFOAM/SubField.H>
32 
33 // * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
34 
35 Foam::scalar Foam::noiseFFT::p0 = 2e-5;
36 
37 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
38 
39 // Construct from pressure field
41 (
42  const scalar deltat,
43  const scalarField& pressure
44 )
45 :
46  scalarField(pressure),
47  deltat_(deltat)
48 {}
49 
50 
51 // Construct from pressure field file name
52 Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
53 :
54  scalarField(),
55  deltat_(0.0)
56 {
57  // Construct control dictionary
58  IFstream pFile(pFileName);
59 
60  // Check pFile stream is OK
61  if (!pFile.good())
62  {
64  (
65  "noiseFFT::noiseFFT(const fileName& pFileName, const label skip)"
66  ) << "Cannot read file " << pFileName
67  << exit(FatalError);
68  }
69 
70  if (skip)
71  {
72  scalar dummyt, dummyp;
73 
74  for (label i=0; i<skip; i++)
75  {
76  pFile >> dummyt;
77 
78  if (!pFile.good() || pFile.eof())
79  {
81  (
82  "noiseFFT::noiseFFT(const fileName& pFileName, "
83  "const label skip)"
84  ) << "Number of points in file " << pFileName
85  << " is less than the number to be skipped = " << skip
86  << exit(FatalError);
87  }
88 
89  pFile >> dummyp;
90  }
91  }
92 
93  scalar t = 0, T = 0;
94  DynamicList<scalar> pData(100000);
95  label i = 0;
96 
97  while (!(pFile >> t).eof())
98  {
99  T = t;
100  pFile >> pData(i++);
101  }
102 
103  deltat_ = T/pData.size();
104 
106 }
107 
108 
109 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
110 
112 {
113  scalarField t(size());
114  forAll (t, i)
115  {
116  t[i] = i*deltat_;
117  }
118 
119  return graph
120  (
121  "p(t)",
122  "t [s]",
123  "p(t) [Pa]",
124  t,
125  *this
126  );
127 }
128 
129 
131 (
132  const label N,
133  const label ni
134 ) const
135 {
136  label windowOffset = N;
137 
138  if ((N + ni*windowOffset) > size())
139  {
140  FatalErrorIn("noiseFFT::window(const label N, const label n) const")
141  << "Requested window is outside set of data" << endl
142  << "number of data = " << size() << endl
143  << "size of window = " << N << endl
144  << "window = " << ni
145  << exit(FatalError);
146  }
147 
148  tmp<scalarField> tpw(new scalarField(N));
149  scalarField& pw = tpw();
150 
151  label offset = ni*windowOffset;
152 
153  forAll (pw, i)
154  {
155  pw[i] = operator[](i + offset);
156  }
157 
158  return tpw;
159 }
160 
161 
163 {
164  scalarField t(N);
165  forAll (t, i)
166  {
167  t[i] = i*deltat_;
168  }
169 
170  scalar T = N*deltat_;
171 
172  return 2*(0.5 - 0.5*cos(2*mathematicalConstant::pi*t/T));
173 }
174 
175 
177 (
178  const tmp<scalarField>& tpn
179 ) const
180 {
181  tmp<scalarField> tPn2
182  (
183  mag
184  (
186  (
187  ReComplexField(tpn),
188  labelList(1, tpn().size())
189  )
190  )
191  );
192 
193  tpn.clear();
194 
195  tmp<scalarField> tPn
196  (
197  new scalarField
198  (
199  scalarField::subField(tPn2(), tPn2().size()/2)
200  )
201  );
202  scalarField& Pn = tPn();
203 
204  Pn *= 2.0/sqrt(scalar(tPn2().size()));
205  Pn[0] /= 2.0;
206 
207  return tPn;
208 }
209 
210 
212 (
213  const label N,
214  const label nw
215 ) const
216 {
217  if (N > size())
218  {
219  FatalErrorIn("noiseFFT::meanPf(const label N, const label nw) const")
220  << "Requested window is outside set of data" << endl
221  << "number of data = " << size() << endl
222  << "size of window = " << N << endl
223  << "window = " << nw
224  << exit(FatalError);
225  }
226 
227  scalarField MeanPf(N/2, 0.0);
228 
229  scalarField Hwf = Hanning(N);
230 
231  for (label wi=0; wi<nw; wi++)
232  {
233  MeanPf += Pf(Hwf*window(N, wi));
234  }
235 
236  MeanPf /= nw;
237 
238  scalarField f(MeanPf.size());
239  scalar deltaf = 1.0/(N*deltat_);
240 
241  forAll (f, i)
242  {
243  f[i] = i*deltaf;
244  }
245 
246  return graph
247  (
248  "P(f)",
249  "f [Hz]",
250  "P(f) [Pa]",
251  f,
252  MeanPf
253  );
254 }
255 
256 
258 (
259  const label N,
260  const label nw
261 ) const
262 {
263  if (N > size())
264  {
265  FatalErrorIn("noiseFFT::RMSmeanPf(const label N, const label nw) const")
266  << "Requested window is outside set of data" << endl
267  << "number of data = " << size() << endl
268  << "size of window = " << N << endl
269  << "window = " << nw
270  << exit(FatalError);
271  }
272 
273  scalarField RMSMeanPf(N/2, 0.0);
274 
275  scalarField Hwf = Hanning(N);
276 
277  for (label wi=0; wi<nw; wi++)
278  {
279  RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
280  }
281 
282  RMSMeanPf = sqrt(RMSMeanPf/nw);
283 
284  scalarField f(RMSMeanPf.size());
285  scalar deltaf = 1.0/(N*deltat_);
286 
287  forAll (f, i)
288  {
289  f[i] = i*deltaf;
290  }
291 
292  return graph
293  (
294  "P(f)",
295  "f [Hz]",
296  "P(f) [Pa]",
297  f,
298  RMSMeanPf
299  );
300 }
301 
302 
304 {
305  return graph
306  (
307  "L(f)",
308  "f [Hz]",
309  "L(f) [dB]",
310  gPf.x(),
311  20*log10(gPf.y()/p0)
312  );
313 }
314 
315 
317 (
318  const graph& gLf,
319  const scalar f1,
320  const scalar fU
321 ) const
322 {
323  const scalarField& f = gLf.x();
324  const scalarField& Lf = gLf.y();
325 
326  scalarField ldelta(Lf.size(), 0.0);
327  scalarField fm(ldelta.size());
328 
329  scalar fratio = cbrt(2.0);
330  scalar deltaf = 1.0/(2*Lf.size()*deltat_);
331 
332  scalar fl = f1/sqrt(fratio);
333  scalar fu = fratio*fl;
334 
335  label istart = label(fl/deltaf);
336  label j = 0;
337 
338  for (label i = istart; i<Lf.size(); i++)
339  {
340  scalar fmi = sqrt(fu*fl);
341 
342  if (fmi > fU + 1) break;
343 
344  if (f[i] >= fu)
345  {
346  fm[j] = fmi;
347  ldelta[j] = 10*log10(ldelta[j]);
348 
349  j++;
350 
351  fl = fu;
352  fu *= fratio;
353  }
354 
355  ldelta[j] += pow(10, Lf[i]/10.0);
356  }
357 
358  fm.setSize(j);
359  ldelta.setSize(j);
360 
361  return graph
362  (
363  "Ldelta",
364  "fm [Hz]",
365  "Ldelta [dB]",
366  fm,
367  ldelta
368  );
369 }
370 
371 
373 (
374  const graph& gPf,
375  const scalar f1,
376  const scalar fU
377 ) const
378 {
379  const scalarField& f = gPf.x();
380  const scalarField& Pf = gPf.y();
381 
382  scalarField pdelta(Pf.size(), 0.0);
383  scalarField fm(pdelta.size());
384 
385  scalar fratio = cbrt(2.0);
386  scalar deltaf = 1.0/(2*Pf.size()*deltat_);
387 
388  scalar fl = f1/sqrt(fratio);
389  scalar fu = fratio*fl;
390 
391  label istart = label(fl/deltaf + 1.0 - SMALL);
392  label j = 0;
393 
394  for (label i = istart; i<Pf.size(); i++)
395  {
396  scalar fmi = sqrt(fu*fl);
397 
398  if (fmi > fU + 1) break;
399 
400  if (f[i] >= fu)
401  {
402  fm[j] = fmi;
403  pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
404 
405  j++;
406 
407  fl = fu;
408  fu *= fratio;
409  }
410 
411  pdelta[j] += sqr(Pf[i]);
412  }
413 
414  fm.setSize(j);
415  pdelta.setSize(j);
416 
417  return graph
418  (
419  "Pdelta",
420  "fm [Hz]",
421  "Pdelta [dB]",
422  fm,
423  pdelta
424  );
425 }
426 
427 
428 Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
429 {
430  const scalarField& Lf = gLf.y();
431 
432  scalar lsum = 0.0;
433 
434  forAll(Lf, i)
435  {
436  lsum += pow(10, Lf[i]/10.0);
437  }
438 
439  lsum = 10*log10(lsum);
440 
441  return lsum;
442 }
443 
444 
445 Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
446 {
447  return p0*pow(10.0, db/20.0);
448 }
449 
450 
452 (
453  const tmp<scalarField>& db
454 ) const
455 {
456  return p0*pow(10.0, db/20.0);
457 }
458 
459 
460 // ************************ vim: set sw=4 sts=4 et: ************************ //