FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

noiseFFT.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "noiseFFT.H"
00027 #include <OpenFOAM/IFstream.H>
00028 #include <OpenFOAM/DynamicList.H>
00029 #include <randomProcesses/fft.H>
00030 #include <OpenFOAM/mathematicalConstants.H>
00031 #include <OpenFOAM/SubField.H>
00032 
00033 // * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * //
00034 
00035 Foam::scalar Foam::noiseFFT::p0 = 2e-5;
00036 
00037 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00038 
00039 // Construct from pressure field
00040 Foam::noiseFFT::noiseFFT
00041 (
00042     const scalar deltat,
00043     const scalarField& pressure
00044 )
00045 :
00046     scalarField(pressure),
00047     deltat_(deltat)
00048 {}
00049 
00050 
00051 // Construct from pressure field file name
00052 Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
00053 :
00054     scalarField(),
00055     deltat_(0.0)
00056 {
00057     // Construct control dictionary
00058     IFstream pFile(pFileName);
00059 
00060     // Check pFile stream is OK
00061     if (!pFile.good())
00062     {
00063         FatalErrorIn
00064         (
00065             "noiseFFT::noiseFFT(const fileName& pFileName, const label skip)"
00066         )   << "Cannot read file " << pFileName
00067             << exit(FatalError);
00068     }
00069 
00070     if (skip)
00071     {
00072         scalar dummyt, dummyp;
00073 
00074         for (label i=0; i<skip; i++)
00075         {
00076             pFile >> dummyt;
00077 
00078             if (!pFile.good() || pFile.eof())
00079             {
00080                 FatalErrorIn
00081                 (
00082                     "noiseFFT::noiseFFT(const fileName& pFileName, "
00083                     "const label skip)"
00084                 )   << "Number of points in file " << pFileName
00085                     << " is less than the number to be skipped = " << skip
00086                     << exit(FatalError);
00087             }
00088 
00089             pFile >> dummyp;
00090         }
00091     }
00092 
00093     scalar t = 0, T = 0;
00094     DynamicList<scalar> pData(100000);
00095     label i = 0;
00096 
00097     while (!(pFile >> t).eof())
00098     {
00099         T = t;
00100         pFile >> pData(i++);
00101     }
00102 
00103     deltat_ = T/pData.size();
00104 
00105     scalarField::operator=(pData.shrink());
00106 }
00107 
00108 
00109 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00110 
00111 Foam::graph Foam::noiseFFT::pt() const
00112 {
00113     scalarField t(size());
00114     forAll (t, i)
00115     {
00116         t[i] = i*deltat_;
00117     }
00118 
00119     return graph
00120     (
00121         "p(t)",
00122         "t [s]",
00123         "p(t) [Pa]",
00124         t,
00125         *this
00126     );
00127 }
00128 
00129 
00130 Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
00131 (
00132     const label N,
00133     const label ni
00134 ) const
00135 {
00136     label windowOffset = N;
00137 
00138     if ((N + ni*windowOffset) > size())
00139     {
00140         FatalErrorIn("noiseFFT::window(const label N, const label n) const")
00141             << "Requested window is outside set of data" << endl
00142             << "number of data = " << size() << endl
00143             << "size of window = " << N << endl
00144             << "window  = " << ni
00145             << exit(FatalError);
00146     }
00147 
00148     tmp<scalarField> tpw(new scalarField(N));
00149     scalarField& pw = tpw();
00150 
00151     label offset = ni*windowOffset;
00152 
00153     forAll (pw, i)
00154     {
00155         pw[i] = operator[](i + offset);
00156     }
00157 
00158     return tpw;
00159 }
00160 
00161 
00162 Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
00163 {
00164     scalarField t(N);
00165     forAll (t, i)
00166     {
00167         t[i] = i*deltat_;
00168     }
00169 
00170     scalar T = N*deltat_;
00171 
00172     return 2*(0.5 - 0.5*cos(2*mathematicalConstant::pi*t/T));
00173 }
00174 
00175 
00176 Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
00177 (
00178     const tmp<scalarField>& tpn
00179 ) const
00180 {
00181     tmp<scalarField> tPn2
00182     (
00183         mag
00184         (
00185             fft::reverseTransform
00186             (
00187                 ReComplexField(tpn),
00188                 labelList(1, tpn().size())
00189             )
00190         )
00191     );
00192 
00193     tpn.clear();
00194 
00195     tmp<scalarField> tPn
00196     (
00197         new scalarField
00198         (
00199             scalarField::subField(tPn2(), tPn2().size()/2)
00200         )
00201     );
00202     scalarField& Pn = tPn();
00203 
00204     Pn *= 2.0/sqrt(scalar(tPn2().size()));
00205     Pn[0] /= 2.0;
00206 
00207     return tPn;
00208 }
00209 
00210 
00211 Foam::graph Foam::noiseFFT::meanPf
00212 (
00213     const label N,
00214     const label nw
00215 ) const
00216 {
00217     if (N > size())
00218     {
00219         FatalErrorIn("noiseFFT::meanPf(const label N, const label nw) const")
00220             << "Requested window is outside set of data" << endl
00221             << "number of data = " << size() << endl
00222             << "size of window = " << N << endl
00223             << "window  = " << nw
00224             << exit(FatalError);
00225     }
00226 
00227     scalarField MeanPf(N/2, 0.0);
00228 
00229     scalarField Hwf = Hanning(N);
00230 
00231     for (label wi=0; wi<nw; wi++)
00232     {
00233         MeanPf += Pf(Hwf*window(N, wi));
00234     }
00235 
00236     MeanPf /= nw;
00237 
00238     scalarField f(MeanPf.size());
00239     scalar deltaf = 1.0/(N*deltat_);
00240 
00241     forAll (f, i)
00242     {
00243         f[i] = i*deltaf;
00244     }
00245 
00246     return graph
00247     (
00248         "P(f)",
00249         "f [Hz]",
00250         "P(f) [Pa]",
00251         f,
00252         MeanPf
00253     );
00254 }
00255 
00256 
00257 Foam::graph Foam::noiseFFT::RMSmeanPf
00258 (
00259     const label N,
00260     const label nw
00261 ) const
00262 {
00263     if (N > size())
00264     {
00265         FatalErrorIn("noiseFFT::RMSmeanPf(const label N, const label nw) const")
00266             << "Requested window is outside set of data" << endl
00267             << "number of data = " << size() << endl
00268             << "size of window = " << N << endl
00269             << "window  = " << nw
00270             << exit(FatalError);
00271     }
00272 
00273     scalarField RMSMeanPf(N/2, 0.0);
00274 
00275     scalarField Hwf = Hanning(N);
00276 
00277     for (label wi=0; wi<nw; wi++)
00278     {
00279         RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
00280     }
00281 
00282     RMSMeanPf = sqrt(RMSMeanPf/nw);
00283 
00284     scalarField f(RMSMeanPf.size());
00285     scalar deltaf = 1.0/(N*deltat_);
00286 
00287     forAll (f, i)
00288     {
00289         f[i] = i*deltaf;
00290     }
00291 
00292     return graph
00293     (
00294         "P(f)",
00295         "f [Hz]",
00296         "P(f) [Pa]",
00297         f,
00298         RMSMeanPf
00299     );
00300 }
00301 
00302 
00303 Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
00304 {
00305     return graph
00306     (
00307         "L(f)",
00308         "f [Hz]",
00309         "L(f) [dB]",
00310         gPf.x(),
00311         20*log10(gPf.y()/p0)
00312     );
00313 }
00314 
00315 
00316 Foam::graph Foam::noiseFFT::Ldelta
00317 (
00318     const graph& gLf,
00319     const scalar f1,
00320     const scalar fU
00321 ) const
00322 {
00323     const scalarField& f = gLf.x();
00324     const scalarField& Lf = gLf.y();
00325 
00326     scalarField ldelta(Lf.size(), 0.0);
00327     scalarField fm(ldelta.size());
00328 
00329     scalar fratio = cbrt(2.0);
00330     scalar deltaf = 1.0/(2*Lf.size()*deltat_);
00331 
00332     scalar fl = f1/sqrt(fratio);
00333     scalar fu = fratio*fl;
00334 
00335     label istart = label(fl/deltaf);
00336     label j = 0;
00337 
00338     for (label i = istart; i<Lf.size(); i++)
00339     {
00340         scalar fmi = sqrt(fu*fl);
00341 
00342         if (fmi > fU + 1) break;
00343 
00344         if (f[i] >= fu)
00345         {
00346             fm[j] = fmi;
00347             ldelta[j] = 10*log10(ldelta[j]);
00348 
00349             j++;
00350 
00351             fl = fu;
00352             fu *= fratio;
00353         }
00354 
00355         ldelta[j] += pow(10, Lf[i]/10.0);
00356     }
00357 
00358     fm.setSize(j);
00359     ldelta.setSize(j);
00360 
00361     return graph
00362     (
00363         "Ldelta",
00364         "fm [Hz]",
00365         "Ldelta [dB]",
00366         fm,
00367         ldelta
00368     );
00369 }
00370 
00371 
00372 Foam::graph Foam::noiseFFT::Pdelta
00373 (
00374     const graph& gPf,
00375     const scalar f1,
00376     const scalar fU
00377 ) const
00378 {
00379     const scalarField& f = gPf.x();
00380     const scalarField& Pf = gPf.y();
00381 
00382     scalarField pdelta(Pf.size(), 0.0);
00383     scalarField fm(pdelta.size());
00384 
00385     scalar fratio = cbrt(2.0);
00386     scalar deltaf = 1.0/(2*Pf.size()*deltat_);
00387 
00388     scalar fl = f1/sqrt(fratio);
00389     scalar fu = fratio*fl;
00390 
00391     label istart = label(fl/deltaf + 1.0 - SMALL);
00392     label j = 0;
00393 
00394     for (label i = istart; i<Pf.size(); i++)
00395     {
00396         scalar fmi = sqrt(fu*fl);
00397 
00398         if (fmi > fU + 1) break;
00399 
00400         if (f[i] >= fu)
00401         {
00402             fm[j] = fmi;
00403             pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
00404 
00405             j++;
00406 
00407             fl = fu;
00408             fu *= fratio;
00409         }
00410 
00411         pdelta[j] += sqr(Pf[i]);
00412     }
00413 
00414     fm.setSize(j);
00415     pdelta.setSize(j);
00416 
00417     return graph
00418     (
00419         "Pdelta",
00420         "fm [Hz]",
00421         "Pdelta [dB]",
00422         fm,
00423         pdelta
00424     );
00425 }
00426 
00427 
00428 Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
00429 {
00430     const scalarField& Lf = gLf.y();
00431 
00432     scalar lsum = 0.0;
00433 
00434     forAll(Lf, i)
00435     {
00436         lsum += pow(10, Lf[i]/10.0);
00437     }
00438 
00439     lsum = 10*log10(lsum);
00440 
00441     return lsum;
00442 }
00443 
00444 
00445 Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
00446 {
00447     return p0*pow(10.0, db/20.0);
00448 }
00449 
00450 
00451 Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
00452 (
00453     const tmp<scalarField>& db
00454 ) const
00455 {
00456     return p0*pow(10.0, db/20.0);
00457 }
00458 
00459 
00460 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines