| 1 | #include "fftwserver.h"
|
|---|
| 2 | #include "FFTW/fftw.h"
|
|---|
| 3 | #include "FFTW/rfftw.h"
|
|---|
| 4 |
|
|---|
| 5 |
|
|---|
| 6 | #define MAXND_FFTW 5
|
|---|
| 7 |
|
|---|
| 8 | class FFTWServerPlan {
|
|---|
| 9 | public:
|
|---|
| 10 | FFTWServerPlan(int n, fftw_direction dir, bool fgreal=false);
|
|---|
| 11 | FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal=false);
|
|---|
| 12 | ~FFTWServerPlan();
|
|---|
| 13 | void Recreate(int n);
|
|---|
| 14 | void Recreate(int nd, int * nxyz);
|
|---|
| 15 |
|
|---|
| 16 | static void FillSizes(const BaseArray & in, int * sz);
|
|---|
| 17 |
|
|---|
| 18 | int _n; // Array dimension for 1-d arrays
|
|---|
| 19 | int _nd; // Nb of dimensions for n-d arrays
|
|---|
| 20 | int _nxyz[MAXND_FFTW]; // Array dimensions for n-d arrays
|
|---|
| 21 | fftw_direction _dir;
|
|---|
| 22 |
|
|---|
| 23 | fftw_plan p;
|
|---|
| 24 | rfftw_plan rp;
|
|---|
| 25 | fftwnd_plan pnd;
|
|---|
| 26 | rfftwnd_plan rpnd;
|
|---|
| 27 |
|
|---|
| 28 | };
|
|---|
| 29 |
|
|---|
| 30 | FFTWServerPlan::FFTWServerPlan(int n, fftw_direction dir, bool fgreal)
|
|---|
| 31 | {
|
|---|
| 32 | if (n < 1)
|
|---|
| 33 | throw ParmError("FFTWServerPlan: Array size <= 0 !");
|
|---|
| 34 | p = NULL;
|
|---|
| 35 | rp = NULL;
|
|---|
| 36 | pnd = NULL;
|
|---|
| 37 | rpnd = NULL;
|
|---|
| 38 | _nd = -1;
|
|---|
| 39 | for(int k=0; k<MAXND_FFTW; k++) _nxyz[k] = -10;
|
|---|
| 40 | _n = n;
|
|---|
| 41 | _dir = dir;
|
|---|
| 42 | if (fgreal) rp = rfftw_create_plan(n, dir, FFTW_ESTIMATE);
|
|---|
| 43 | else p = fftw_create_plan(n, dir, FFTW_ESTIMATE);
|
|---|
| 44 | }
|
|---|
| 45 |
|
|---|
| 46 | FFTWServerPlan::FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal)
|
|---|
| 47 | {
|
|---|
| 48 | int k;
|
|---|
| 49 | if (nd > MAXND_FFTW)
|
|---|
| 50 | throw ParmError("FFTWServerPlan: Array rank (nd) > MAXND_FFTW !");
|
|---|
| 51 | p = NULL;
|
|---|
| 52 | rp = NULL;
|
|---|
| 53 | pnd = NULL;
|
|---|
| 54 | rpnd = NULL;
|
|---|
| 55 | _n = -10;
|
|---|
| 56 | _nd = nd;
|
|---|
| 57 | for(k=0; k<nd; k++) {
|
|---|
| 58 | if (nxyz[k] < 1)
|
|---|
| 59 | throw ParmError("FFTWServerPlan: One of the Array size <= 0 !");
|
|---|
| 60 | _nxyz[k] = nxyz[k];
|
|---|
| 61 | }
|
|---|
| 62 | for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
|
|---|
| 63 | _dir = dir;
|
|---|
| 64 | if (fgreal) rpnd = rfftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
|
|---|
| 65 | else pnd = fftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
|
|---|
| 66 | }
|
|---|
| 67 |
|
|---|
| 68 | FFTWServerPlan::~FFTWServerPlan()
|
|---|
| 69 | {
|
|---|
| 70 | if (p) fftw_destroy_plan(p);
|
|---|
| 71 | if (rp) rfftw_destroy_plan(rp);
|
|---|
| 72 | if (pnd) fftwnd_destroy_plan(pnd);
|
|---|
| 73 | if (rpnd) rfftwnd_destroy_plan(rpnd);
|
|---|
| 74 | }
|
|---|
| 75 |
|
|---|
| 76 | void
|
|---|
| 77 | FFTWServerPlan::Recreate(int n)
|
|---|
| 78 | {
|
|---|
| 79 | if (n < 1)
|
|---|
| 80 | throw ParmError("FFTWServerPlan::Recreate(n) n < 0 !");
|
|---|
| 81 | if (_nd > 0)
|
|---|
| 82 | throw ParmError("FFTWServerPlan::Recreate(n) Multi-dimensional plan ! > 0 !");
|
|---|
| 83 | if (n == _n) return;
|
|---|
| 84 | _n = n;
|
|---|
| 85 | if (p) {
|
|---|
| 86 | fftw_destroy_plan(p);
|
|---|
| 87 | p = fftw_create_plan(n, _dir, FFTW_ESTIMATE);
|
|---|
| 88 | }
|
|---|
| 89 | else {
|
|---|
| 90 | rfftw_destroy_plan(rp);
|
|---|
| 91 | rp = rfftw_create_plan(n, _dir, FFTW_ESTIMATE);
|
|---|
| 92 | }
|
|---|
| 93 | }
|
|---|
| 94 |
|
|---|
| 95 | void
|
|---|
| 96 | FFTWServerPlan::Recreate(int nd, int * nxyz)
|
|---|
| 97 | {
|
|---|
| 98 | if (_n > 0)
|
|---|
| 99 | throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) 1-dimensional plan !");
|
|---|
| 100 | int k;
|
|---|
| 101 | if (nd == _nd) {
|
|---|
| 102 | bool samepl = true;
|
|---|
| 103 | for (int k=0; k<nd; k++)
|
|---|
| 104 | if (nxyz[k] != _nxyz[k]) samepl = false;
|
|---|
| 105 | if (samepl) return;
|
|---|
| 106 | }
|
|---|
| 107 | if (nd > MAXND_FFTW)
|
|---|
| 108 | throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) Array rank (nd) > MAXND_FFTW !");
|
|---|
| 109 | _nd = nd;
|
|---|
| 110 | for(k=0; k<nd; k++) {
|
|---|
| 111 | if (nxyz[k] < 1)
|
|---|
| 112 | throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) One of the Array size <= 0 !");
|
|---|
| 113 | _nxyz[k] = nxyz[k];
|
|---|
| 114 | }
|
|---|
| 115 | for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
|
|---|
| 116 | if (pnd) {
|
|---|
| 117 | fftwnd_destroy_plan(pnd);
|
|---|
| 118 | pnd = fftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
|
|---|
| 119 | }
|
|---|
| 120 | else {
|
|---|
| 121 | rfftwnd_destroy_plan(rpnd);
|
|---|
| 122 | rpnd = rfftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
|
|---|
| 123 | }
|
|---|
| 124 |
|
|---|
| 125 | }
|
|---|
| 126 |
|
|---|
| 127 | void
|
|---|
| 128 | FFTWServerPlan::FillSizes(const BaseArray & in, int * sz)
|
|---|
| 129 | {
|
|---|
| 130 | int k1 = 0;
|
|---|
| 131 | int k2 = 0;
|
|---|
| 132 | for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
|---|
| 133 | sz[k2] = in.Size(k1); k2++;
|
|---|
| 134 | }
|
|---|
| 135 | }
|
|---|
| 136 |
|
|---|
| 137 | /* --Methode-- */
|
|---|
| 138 | FFTWServer::FFTWServer()
|
|---|
| 139 | : FFTServerInterface("FFTServer using FFTW package")
|
|---|
| 140 | , ckR8("FFTWServer - ", true, false)
|
|---|
| 141 |
|
|---|
| 142 | {
|
|---|
| 143 | _p1df = NULL;
|
|---|
| 144 | _p1db = NULL;
|
|---|
| 145 | _pndf = NULL;
|
|---|
| 146 | _pndb = NULL;
|
|---|
| 147 |
|
|---|
| 148 | _p1drf = NULL;
|
|---|
| 149 | _p1drb = NULL;
|
|---|
| 150 | _pndrf = NULL;
|
|---|
| 151 | _pndrb = NULL;
|
|---|
| 152 | }
|
|---|
| 153 |
|
|---|
| 154 |
|
|---|
| 155 | /* --Methode-- */
|
|---|
| 156 | FFTWServer::~FFTWServer()
|
|---|
| 157 | {
|
|---|
| 158 | if (_p1df) delete _p1df ;
|
|---|
| 159 | if (_p1db) delete _p1db ;
|
|---|
| 160 | if (_pndf) delete _pndf ;
|
|---|
| 161 | if (_pndb) delete _pndb ;
|
|---|
| 162 |
|
|---|
| 163 | if (_p1drf) delete _p1drf ;
|
|---|
| 164 | if (_p1drb) delete _p1drb ;
|
|---|
| 165 | if (_pndrf) delete _pndrf ;
|
|---|
| 166 | if (_pndrb) delete _pndrb ;
|
|---|
| 167 | }
|
|---|
| 168 |
|
|---|
| 169 | /* --Methode-- */
|
|---|
| 170 | FFTServerInterface * FFTWServer::Clone()
|
|---|
| 171 | {
|
|---|
| 172 | return (new FFTWServer) ;
|
|---|
| 173 | }
|
|---|
| 174 |
|
|---|
| 175 | /* --Methode-- */
|
|---|
| 176 | void
|
|---|
| 177 | FFTWServer::FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
|
|---|
| 178 | {
|
|---|
| 179 | int rank = ckR8.CheckResize(in, out);
|
|---|
| 180 | if (rank == 1) { // One dimensional transform
|
|---|
| 181 | if (_p1df) _p1df->Recreate(in.Size());
|
|---|
| 182 | else _p1df = new FFTWServerPlan(in.Size(), FFTW_FORWARD, false);
|
|---|
| 183 | fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
|---|
| 184 | }
|
|---|
| 185 | else { // Multi dimensional
|
|---|
| 186 | if (in.NbDimensions() > MAXND_FFTW)
|
|---|
| 187 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
|
|---|
| 188 | int sz[MAXND_FFTW];
|
|---|
| 189 | FFTWServerPlan::FillSizes(in, sz);
|
|---|
| 190 | // int k1 = 0;
|
|---|
| 191 | // int k2 = 0;
|
|---|
| 192 | // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
|---|
| 193 | // sz[k2] = in.Size(k1); k2++;
|
|---|
| 194 | // }
|
|---|
| 195 | if (_pndf) _pndf->Recreate(in.NbDimensions(), sz);
|
|---|
| 196 | else _pndf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_FORWARD, false);
|
|---|
| 197 | fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
|---|
| 198 | }
|
|---|
| 199 | if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
|
|---|
| 200 | return;
|
|---|
| 201 | }
|
|---|
| 202 |
|
|---|
| 203 | /* --Methode-- */
|
|---|
| 204 | void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
|
|---|
| 205 | {
|
|---|
| 206 | int rank = ckR8.CheckResize(in, out);
|
|---|
| 207 | if (rank == 1) { // One dimensional transform
|
|---|
| 208 | if (_p1db) _p1db->Recreate(in.Size());
|
|---|
| 209 | else _p1db = new FFTWServerPlan(in.Size(), FFTW_BACKWARD, false);
|
|---|
| 210 | fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
|---|
| 211 | }
|
|---|
| 212 | else { // Multi dimensional
|
|---|
| 213 | if (in.NbDimensions() > MAXND_FFTW)
|
|---|
| 214 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
|
|---|
| 215 | int sz[MAXND_FFTW];
|
|---|
| 216 | FFTWServerPlan::FillSizes(in, sz);
|
|---|
| 217 | // int k1 = 0;
|
|---|
| 218 | // int k2 = 0;
|
|---|
| 219 | // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
|---|
| 220 | // sz[k2] = in.Size(k1); k2++;
|
|---|
| 221 | // }
|
|---|
| 222 | if (_pndb) _pndb->Recreate(in.NbDimensions(), sz);
|
|---|
| 223 | else _pndb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_BACKWARD, false);
|
|---|
| 224 | fftwnd_one(_pndb->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
|---|
| 225 | }
|
|---|
| 226 |
|
|---|
| 227 | return;
|
|---|
| 228 | }
|
|---|
| 229 |
|
|---|
| 230 |
|
|---|
| 231 | void FFTWServer::FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out)
|
|---|
| 232 | {
|
|---|
| 233 | int rank = ckR8.CheckResize(in, out);
|
|---|
| 234 | if (rank == 1) { // One dimensional transform
|
|---|
| 235 | if (_p1drf) _p1drf->Recreate(in.Size());
|
|---|
| 236 | else _p1drf = new FFTWServerPlan(in.Size(), FFTW_REAL_TO_COMPLEX, true);
|
|---|
| 237 | TArray<r_8> outtemp;
|
|---|
| 238 | outtemp.ReSize(in);
|
|---|
| 239 | rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outtemp.Data()));
|
|---|
| 240 | ReShapetoCompl(outtemp, out);
|
|---|
| 241 | }
|
|---|
| 242 | else { // Multi dimensional
|
|---|
| 243 | if (in.NbDimensions() > MAXND_FFTW)
|
|---|
| 244 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
|
|---|
| 245 | int sz[MAXND_FFTW];
|
|---|
| 246 | FFTWServerPlan::FillSizes(in, sz);
|
|---|
| 247 | // int k1 = 0;
|
|---|
| 248 | // int k2 = 0;
|
|---|
| 249 | // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
|---|
| 250 | // sz[k2] = in.Size(k1); k2++;
|
|---|
| 251 | // }
|
|---|
| 252 | if (_pndrf) _pndrf->Recreate(in.NbDimensions(), sz);
|
|---|
| 253 | else _pndrf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_REAL_TO_COMPLEX, true);
|
|---|
| 254 | rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) ,
|
|---|
| 255 | (fftw_complex *)(out.Data()) );
|
|---|
| 256 | }
|
|---|
| 257 | if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
|
|---|
| 258 | return;
|
|---|
| 259 |
|
|---|
| 260 | }
|
|---|
| 261 |
|
|---|
| 262 |
|
|---|
| 263 |
|
|---|
| 264 | void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out,
|
|---|
| 265 | bool usoutsz)
|
|---|
| 266 | {
|
|---|
| 267 |
|
|---|
| 268 | int rank = ckR8.CheckResize(in, out, usoutsz);
|
|---|
| 269 | if (rank == 1) { // One dimensional transform
|
|---|
| 270 | TArray<r_8> intemp;
|
|---|
| 271 | intemp.ReSize(out);
|
|---|
| 272 | if (_p1drb) _p1drb->Recreate(out.Size());
|
|---|
| 273 | else _p1drb = new FFTWServerPlan(out.Size(), FFTW_COMPLEX_TO_REAL, true);
|
|---|
| 274 |
|
|---|
| 275 | ReShapetoReal(in, intemp);
|
|---|
| 276 | // cerr << " DEBUG-FFTWServer::FFTBackward() in = \n" << in << endl;
|
|---|
| 277 | // cerr << " DEBUG-FFTWServer::FFTBackward() intemp = \n" << intemp << endl;
|
|---|
| 278 | rfftw_one(_p1drb->rp, (fftw_real *)(intemp.Data()) , (fftw_real *)(out.Data()));
|
|---|
| 279 |
|
|---|
| 280 | }
|
|---|
| 281 | else { // Multi dimensional
|
|---|
| 282 | if (in.NbDimensions() > MAXND_FFTW)
|
|---|
| 283 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, r_8 ) rank > MAXND_FFTW !");
|
|---|
| 284 | int sz[MAXND_FFTW];
|
|---|
| 285 | FFTWServerPlan::FillSizes(out, sz);
|
|---|
| 286 | // int k1 = 0;
|
|---|
| 287 | // int k2 = 0;
|
|---|
| 288 | // for(k1=out.NbDimensions()-1; k1>=0; k1--) {
|
|---|
| 289 | // sz[k2] = out.Size(k1); k2++;
|
|---|
| 290 | // }
|
|---|
| 291 | if (_pndrb) _pndrb->Recreate(in.NbDimensions(), sz);
|
|---|
| 292 | else _pndrb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_COMPLEX_TO_REAL, true);
|
|---|
| 293 |
|
|---|
| 294 | rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) );
|
|---|
| 295 | }
|
|---|
| 296 | return;
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 |
|
|---|
| 300 |
|
|---|
| 301 | /* --Methode-- */
|
|---|
| 302 | void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa)
|
|---|
| 303 | {
|
|---|
| 304 | TVector< complex<r_8> > in(ina);
|
|---|
| 305 | TVector< r_8> out(outa);
|
|---|
| 306 | int n = in.NElts();
|
|---|
| 307 | r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
|
|---|
| 308 | sa_size_t ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
|
|---|
| 309 | ncs = 2*n-1 : ncs = 2*n-2;
|
|---|
| 310 |
|
|---|
| 311 | if (out.NElts() != ncs) {
|
|---|
| 312 | cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
|
|---|
| 313 | << " out.NElts()= " << out.NElts() << endl;
|
|---|
| 314 | throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
|
|---|
| 315 | }
|
|---|
| 316 | // if (ncs == 2*n-2) {
|
|---|
| 317 | out(0) = in(0).real();
|
|---|
| 318 | for(int k=1; k<n; k++) {
|
|---|
| 319 | out(k) = in(k).real();
|
|---|
| 320 | out(ncs-k) = in(k).imag();
|
|---|
| 321 | }
|
|---|
| 322 | if (ncs == 2*n-2)
|
|---|
| 323 | out(ncs-1) = in(n-1).real();
|
|---|
| 324 |
|
|---|
| 325 | // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
|
|---|
| 326 | }
|
|---|
| 327 |
|
|---|
| 328 |
|
|---|
| 329 | /* --Methode-- */
|
|---|
| 330 | void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
|
|---|
| 331 | {
|
|---|
| 332 | TVector< r_8> in(ina);
|
|---|
| 333 | TVector< complex<r_8> > out(outa);
|
|---|
| 334 |
|
|---|
| 335 | sa_size_t n = in.NElts();
|
|---|
| 336 | sa_size_t ncs = n/2+1;
|
|---|
| 337 | sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
|
|---|
| 338 | if (out.NElts() != ncs)
|
|---|
| 339 | throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
|
|---|
| 340 |
|
|---|
| 341 | out(0) = complex<r_8> (in(0),0.);
|
|---|
| 342 | for(int k=1; k<ncs; k++)
|
|---|
| 343 | out(k) = complex<r_8>(in(k), in(n-k));
|
|---|
| 344 | if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
|
|---|
| 345 |
|
|---|
| 346 | }
|
|---|
| 347 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
