src/randomProcesses/noise/noiseFFT.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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 "mathematicalConstants.H"
  31 #include "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
  40 Foam::noiseFFT::noiseFFT
  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     {
  63         FatalErrorIn
  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             {
  80                 FatalErrorIn
  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
 105     scalarField::operator=(pData.shrink());
 106 }
 107
 108
 109 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 110
 111 Foam::graph Foam::noiseFFT::pt() const
 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
 130 Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
 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
 162 Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
 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
 176 Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
 177 (
 178     const tmp<scalarField>& tpn
 179 ) const
 180 {
 181     tmp<scalarField> tPn2
 182     (
 183         mag
 184         (
 185             fft::reverseTransform
 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
 211 Foam::graph Foam::noiseFFT::meanPf
 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
 257 Foam::graph Foam::noiseFFT::RMSmeanPf
 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
 303 Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
 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
 316 Foam::graph Foam::noiseFFT::Ldelta
 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
 372 Foam::graph Foam::noiseFFT::Pdelta
 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
 451 Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
 452 (
 453     const tmp<scalarField>& db
 454 ) const
 455 {
 456     return p0*pow(10.0, db/20.0);
 457 }
 458
 459
 460 // ************************************************************************* //