40 void ProSHADE_internal_sphericalHarmonics::allocateComputationMemory ( proshade_unsign band, proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& outputReal, proshade_double*& outputImag, proshade_double*& shWeights,
double**& tableSpace, proshade_double*& tableSpaceHelper, fftw_complex*& workspace )
43 proshade_unsign oneDimmension = 2 * band;
46 inputReal =
new proshade_double [oneDimmension * oneDimmension];
47 inputImag =
new proshade_double [oneDimmension * oneDimmension];
50 outputReal =
new proshade_double [oneDimmension * oneDimmension];
51 outputImag =
new proshade_double [oneDimmension * oneDimmension];
54 shWeights =
new proshade_double [band * 4];
55 workspace =
new fftw_complex [( 8 * band * band ) + ( 10 * band )];
58 tableSpaceHelper =
new proshade_double [
static_cast<proshade_unsign
> ( Reduced_Naive_TableSize ( band, band ) +
59 Reduced_SpharmonicTableSize ( band, band ) )];
90 rres =
reinterpret_cast<proshade_double*
> ( workspace );
91 ires = rres + ( oDim * oDim );
92 fltres = ires + ( oDim * oDim );
93 scratchpad = fltres + ( oDim / 2 );
118 fftw_iodim howmany_dims[1];
121 int howmany_rank = 1;
123 dims[0].n =
static_cast<int> ( band * 2 );
125 dims[0].os =
static_cast<int> ( band * 2 );
127 howmany_dims[0].n =
static_cast<int> ( band * 2 );
128 howmany_dims[0].is =
static_cast<int> ( band * 2 );
129 howmany_dims[0].os = 1;
132 fftPlan = fftw_plan_guru_split_dft ( rank,
143 dctPlan = fftw_plan_r2r_1d (
static_cast<int> ( band * 2 ),
145 scratchpad +
static_cast<int> ( band * 2 ),
171 void ProSHADE_internal_sphericalHarmonics::releaseSphericalMemory ( proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& outputReal, proshade_double*& outputImag,
double*& tableSpaceHelper,
double**& tableSpace,
double*& shWeights, fftw_complex*& workspace, fftw_plan& fftPlan, fftw_plan& dctPlan, proshade_unsign band )
178 delete[] tableSpaceHelper;
183 tableSpaceHelper = NULL;
189 fftw_destroy_plan ( dctPlan );
190 fftw_destroy_plan ( fftPlan );
218 void ProSHADE_internal_sphericalHarmonics::initialiseAllMemory ( proshade_unsign band, proshade_double*& inputReal, proshade_double*& inputImag, proshade_double*& outputReal, proshade_double*& outputImag,
double*& shWeights,
double**& tableSpace,
double*& tableSpaceHelper, fftw_complex*& workspace, proshade_double*& rres, proshade_double*& ires, proshade_double*& fltres, proshade_double*& scratchpad, fftw_plan& fftPlan, fftw_plan& dctPlan )
221 proshade_unsign oneDim = band * 2;
224 allocateComputationMemory ( band, inputReal, inputImag, outputReal, outputImag, shWeights, tableSpace, tableSpaceHelper, workspace );
230 tableSpace = SemiNaive_Naive_Pml_Table ( band, band, tableSpaceHelper,
reinterpret_cast<double*
> ( workspace ) );
233 makeweights ( band, shWeights );
236 initialiseFFTWPlans ( band, fftPlan, dctPlan, inputReal, inputImag, rres, ires, scratchpad );
261 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( oneDim * oneDim ); iter++ )
263 inputReal[iter] = mappedData[iter];
264 inputImag[iter] = 0.0;
268 fftw_execute_split_dft ( fftPlan, inputReal, inputImag, rres, ires ) ;
271 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( oneDim * oneDim ); iter++ )
273 rres[iter] *= normCoeff;
274 ires[iter] *= normCoeff;
301 void ProSHADE_internal_sphericalHarmonics::computeSphericalTransformCoeffs ( proshade_unsign band, proshade_double*& rdataptr, proshade_double*& idataptr, proshade_double*& outputReal, proshade_double*& outputImag, proshade_double*& rres, proshade_double*& ires, proshade_double*& fltres, proshade_double*& scratchpad,
double**& tablePml,
double*& shWeights, fftw_plan& dctPlan )
304 rdataptr = outputReal;
305 idataptr = outputImag;
306 for ( proshade_unsign bandIter = 0; bandIter < band; bandIter++ )
309 SemiNaiveReduced ( rres + ( bandIter * ( band * 2 ) ),
319 memcpy ( rdataptr, fltres,
sizeof(proshade_double) * ( band - bandIter ) );
320 rdataptr += band - bandIter;
323 SemiNaiveReduced ( ires + ( bandIter * ( band * 2 ) ),
333 memcpy ( idataptr, fltres,
sizeof(proshade_double) * ( band - bandIter ) );
334 idataptr += band - bandIter;
356 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( (band * 2) * (band * 2) ); iter++ )
358 shArray[iter][0] = outputReal[iter];
359 shArray[iter][1] = outputImag[iter];
363 proshade_double powerOne = 1.0;
364 proshade_unsign hlp1 = 0;
365 proshade_unsign hlp2 = 0;
366 for ( proshade_signed order = 1; order < static_cast<proshade_signed> ( band ); order++)
369 for ( proshade_signed bandIter = order; bandIter < static_cast<proshade_signed> ( band ); bandIter++)
371 hlp1 = seanindex ( order, bandIter, band );
372 hlp2 = seanindex ( -order, bandIter, band );
374 shArray[hlp2][0] = powerOne *
static_cast<proshade_double
> ( outputReal[hlp1] );
375 shArray[hlp2][1] = -powerOne *
static_cast<proshade_double
> ( outputImag[hlp1] );
398 proshade_double *inputReal = NULL, *inputImag = NULL, *outputReal = NULL, *outputImag = NULL;
399 double *shWeights = NULL, *tableSpaceHelper = NULL;
400 double** tablePml = NULL;
401 fftw_complex* workspace = NULL;
402 proshade_unsign oneDim =
static_cast<proshade_unsign
> ( band * 2 );
403 proshade_double normCoeff = ( 1.0 / (
static_cast<proshade_double
> ( band * 2 ) ) ) * sqrt( 2.0 * M_PI );
406 for ( proshade_unsign i = 0; i < ( 2 * band * 2 * band); i++ )
413 proshade_double *rres = NULL, *ires = NULL, *fltres = NULL, *scratchpad = NULL, *rdataptr = NULL, *idataptr = NULL;
416 fftw_plan fftPlan = NULL;
417 fftw_plan dctPlan = NULL;
420 initialiseAllMemory ( band, inputReal, inputImag, outputReal, outputImag, shWeights, tablePml, tableSpaceHelper, workspace,
421 rres, ires, fltres, scratchpad, fftPlan, dctPlan );
427 computeSphericalTransformCoeffs ( band, rdataptr, idataptr, outputReal, outputImag, rres, ires, fltres, scratchpad, tablePml, shWeights, dctPlan );
433 releaseSphericalMemory ( inputReal, inputImag, outputReal, outputImag, tableSpaceHelper, tablePml, shWeights, workspace, fftPlan, dctPlan, band );