Source
308
308
if(x%2 != 0 || y%2 !=0)
309
309
throw(AipsError("Programmer error: FFT2D does not deal with odd numbers on x-y plane"));
310
310
fftShift(out, x, y, true);
311
311
doFFT(out, x, y, toFreq);
312
312
fftShift(out, x, y, toFreq);
313
313
314
314
}
315
315
void FFT2D::doFFT(DComplex*& out, Long x, Long y, Bool toFreq){
316
316
if(useFFTW_p){
317
317
//Will need to seperate the plan from the execute if we want to run this in multiple threads
318
-
Int dim[2]={Int(x), Int(y)};
318
+
Int dim[2]={Int(y), Int(x)};
319
319
if(toFreq){
320
320
321
321
planC2CD_p=fftw_plan_dft(2, dim, reinterpret_cast<fftw_complex *>(out), reinterpret_cast<fftw_complex *>(out), FFTW_FORWARD, FFTW_ESTIMATE);
322
322
323
323
//fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
324
324
}
325
325
else{
326
326
planC2CD_p=fftw_plan_dft(2, dim, reinterpret_cast<fftw_complex *>(out), reinterpret_cast<fftw_complex *>(out), FFTW_BACKWARD, FFTW_ESTIMATE);
327
327
// fft1_p.plan_c2c_backward(IPosition(2, x, y), out);
328
328
}
329
329
fftw_execute(planC2CD_p);
330
330
331
331
}
332
332
else{
333
333
throw(AipsError("Double precision FFT with FFTPack is not implemented"));
334
334
}
335
335
}
336
336
void FFT2D::doFFT(Complex*& out, Long x, Long y, Bool toFreq){
337
337
if(useFFTW_p){
338
338
//Will need to seperate the plan from the execute if we want to run this in multiple threads
339
-
Int dim[2]={Int(x), Int(y)};
339
+
Int dim[2]={Int(y), Int(x)};
340
340
if(toFreq){
341
341
342
342
planC2C_p=fftwf_plan_dft(2, dim, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out), FFTW_FORWARD, FFTW_ESTIMATE);
343
343
344
344
//fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
345
345
}
346
346
else{
347
347
planC2C_p=fftwf_plan_dft(2, dim, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out), FFTW_BACKWARD, FFTW_ESTIMATE);
348
348
// fft1_p.plan_c2c_backward(IPosition(2, x, y), out);
349
349
}
358
358
wsave_p.resize(2*x1*y1+15);
359
359
lsav_p=2*x1*y1+15;
360
360
Float *wsaveptr=wsave_p.data();
361
361
FFTPack::cfft2i(x1, y1, wsaveptr, lsav_p, ier);
362
362
}
363
363
std::vector<Float> work(2*x1*y1);
364
364
Int lenwrk=2*x1*y1;
365
365
Float* workptr=work.data();
366
366
Float* wsaveptr=wsave_p.data();
367
367
if(toFreq)
368
-
FFTPack::cfft2f(x1, y1, x1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
368
+
FFTPack::cfft2f(x1, x1, y1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
369
369
else
370
-
FFTPack::cfft2b(x1, y1, x1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
370
+
FFTPack::cfft2b(x1, x1, y1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
371
371
}
372
372
}
373
373
void FFT2D::doFFT(Complex*& out, Float*& in, Long x, Long y){
374
374
if(useFFTW_p){
375
-
Int dim[2]={Int(x), Int(y)};
375
+
Int dim[2]={Int(y), Int(x)};
376
376
377
377
planR2C_p=fftwf_plan_dft_r2c(2, dim, in, reinterpret_cast<fftwf_complex *>(out), FFTW_ESTIMATE);
378
378
379
379
//fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
380
380
381
381
fftwf_execute(planR2C_p);
382
382
383
383
}
384
384
else{
385
385
/*