00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
00034
00035 Foam::scalar Foam::noiseFFT::p0 = 2e-5;
00036
00037
00038
00039
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
00052 Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
00053 :
00054 scalarField(),
00055 deltat_(0.0)
00056 {
00057
00058 IFstream pFile(pFileName);
00059
00060
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
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