ENH: FFT - updated for fft of real data

This commit is contained in:
Andrew Heather
2017-11-16 18:05:56 +00:00
parent 929ed3ca02
commit b2f002ba45
4 changed files with 68 additions and 29 deletions

View File

@ -34,6 +34,56 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field)
{
// Copy of input field for use by fftw
// - fftw requires non-const access to input and output...
scalarField fftInOut(field);
const label n = fftInOut.size();
const label nBy2 = n/2;
// Using real to half-complex fftw 'kind'
fftw_plan plan = fftw_plan_r2r_1d
(
n,
fftInOut.data(),
fftInOut.data(),
FFTW_R2HC,
FFTW_ESTIMATE
);
fftw_execute(plan);
fftw_destroy_plan(plan);
// field[0] = DC component
tmp<complexField> tresult(new complexField(nBy2 + 1));
complexField& result = tresult.ref();
result[0].Re() = fftInOut[0];
result[nBy2].Re() = fftInOut[nBy2];
for (label i = 1; i < nBy2; ++i)
{
result[i].Re() = fftInOut[i];
result[i].Im() = fftInOut[n - i];
}
return tresult;
}
Foam::tmp<Foam::complexField> fft::realTransform1D
(
const tmp<scalarField>& tfield
)
{
tmp<complexField> tresult = realTransform1D(tfield());
tfield.clear();
return tresult;
}
void fft::transform
(
complexField& field,

View File

@ -44,7 +44,6 @@ SourceFiles
#define fft_H
#include "complexFields.H"
#include "UList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,6 +62,19 @@ public:
};
//- Transform real-value data
// - uses the fftw real to half-complex method
// - result size is field.size()/2 + 1
static tmp<complexField> realTransform1D(const scalarField& field);
//- Transform real-value data
// - uses the fftw real to half-complex method
// - result size is field.size()/2 + 1
static tmp<complexField> realTransform1D(const tmp<scalarField>& field);
//- Transform compex-value data
static void transform
(
complexField& field,

View File

@ -26,7 +26,6 @@ License
#include "noiseFFT.H"
#include "IFstream.H"
#include "DynamicList.H"
#include "SubField.H"
#include "mathematicalConstants.H"
#include "HashSet.H"
#include "fft.H"
@ -240,33 +239,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
const tmp<scalarField>& tpn
) const
{
// Calculate the 2-sided fft
// Note: result not scaled
tmp<scalarField> tPn2
(
mag
(
fft::reverseTransform
(
ReComplexField(tpn),
List<int>(1, tpn().size())
)
)
);
tpn.clear();
// Only storing the positive half
// Note: storing (N/2) values, DC component at position (0)
tmp<scalarField> tPn
(
new scalarField
(
scalarField::subField(tPn2(), tPn2().size()/2 + 1)
)
);
return tPn;
return mag(fft::realTransform1D(tpn));
}

View File

@ -23,6 +23,10 @@ License
\*---------------------------------------------------------------------------*/
#include "SubField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::windowModel::apply
(