src/randomProcesses/noise/noiseFFT.C

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