Line data Source code
1 : //# FFT2D.cc: implementation of FFT2D
2 : //# Copyright (C) 2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $kgolap$
28 : //DEDICATED TO HONGLIN YE
29 :
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/Array.h>
32 : #include <casacore/casa/OS/HostInfo.h>
33 : #include <synthesis/Utilities/FFT2D.h>
34 : #include <casacore/lattices/Lattices/Lattice.h>
35 : #ifdef _OPENMP
36 : #include <omp.h>
37 : #endif
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 : // Utility function to try to give as much info as possible - CAS-12604
42 0 : void throw_programmer_error(Long nx_p, Long ny_p, Long x, Long y, const char *file,
43 : int line)
44 : {
45 0 : AipsError exc;
46 0 : ostringstream msg;
47 0 : msg << "Programmer error: called FFT2D with wrong size. In " << file << ":" << line
48 0 : << ", with nx_p: " << nx_p << ", ny_p: " << ny_p << ", x: " << x << ", y: " << y
49 0 : << "\n Stack trace: " << exc.getStackTrace();
50 0 : exc.setMessage(msg.str());
51 0 : throw(exc);
52 0 : }
53 :
54 118 : FFT2D::FFT2D(Bool useFFTW):planC2C_forw_p(nullptr),planC2C_back_p(nullptr), planR2C_p(nullptr), planC2CD_forw_p(nullptr), planC2CD_back_p(nullptr), useFFTW_p(useFFTW), wsave_p(0), lsav_p(0), nx_p(-1), ny_p(-1) {
55 118 : if(useFFTW_p){
56 118 : Int numThreads=HostInfo::numCPUs(true);
57 : #ifdef _OPENMP
58 118 : numThreads=omp_get_max_threads();
59 : #endif
60 118 : fftwf_init_threads();
61 118 : fftwf_plan_with_nthreads(numThreads);
62 : ///For double precision
63 118 : fftw_init_threads();
64 118 : fftw_plan_with_nthreads(numThreads);
65 : }
66 :
67 118 : }
68 114 : FFT2D::~FFT2D(){
69 114 : if(useFFTW_p){
70 114 : if(planC2CD_forw_p)
71 0 : fftw_destroy_plan(planC2CD_forw_p);
72 114 : if(planC2C_forw_p)
73 0 : fftwf_destroy_plan(planC2C_forw_p);
74 114 : if(planR2C_p)
75 0 : fftwf_destroy_plan(planR2C_p);
76 114 : if(planC2CD_back_p)
77 0 : fftw_destroy_plan(planC2CD_back_p);
78 114 : if(planC2C_back_p)
79 0 : fftwf_destroy_plan(planC2C_back_p);
80 : ////Have to leak the cleanup part as it is thread unsafe to perform this, see CAS-12486
81 : // fftw_cleanup_threads();
82 : // fftw_cleanup();
83 : // fftwf_cleanup_threads();
84 : // fftwf_cleanup();
85 :
86 114 : planC2CD_forw_p=nullptr;
87 114 : planC2C_forw_p=nullptr;
88 114 : planC2CD_back_p=nullptr;
89 114 : planC2C_back_p=nullptr;
90 : }
91 114 : }
92 :
93 3 : FFT2D& FFT2D::operator=(const FFT2D& other){
94 3 : if(this != &other){
95 : /*planC2C_forw_p=other.planC2C_forw_p;
96 : planR2C_p=other.planR2C_p;
97 : planC2CD_forw_p=other.planC2CD_forw_p;
98 : planC2C_back_p=other.planC2C_back_p;
99 : planC2CD_back_p=other.planC2CD_back_p;
100 : */
101 3 : planC2C_forw_p=nullptr;
102 3 : planR2C_p=nullptr;
103 3 : planC2CD_forw_p=nullptr;
104 3 : planC2C_back_p=nullptr;
105 3 : planC2CD_back_p=nullptr;
106 : //cerr << "copy: planc2cd " << planC2CD_back_p << endl;
107 3 : useFFTW_p=other.useFFTW_p;
108 3 : wsave_p.resize(other.wsave_p.size());
109 3 : wsave_p=other.wsave_p;
110 3 : lsav_p=other.lsav_p;
111 3 : nx_p=-1;
112 3 : ny_p=-1;
113 :
114 : }
115 3 : return *this;
116 : }
117 :
118 0 : void FFT2D::r2cFFT(Lattice<Complex>& out, Lattice<Float>& in){
119 :
120 0 : IPosition shp=in.shape();
121 0 : if(shp.nelements() <2)
122 0 : throw(AipsError("Lattice has to be 2 dimensional to use FFT2D"));
123 0 : Long x= in.shape()(0);
124 0 : Long y=in.shape()(1);
125 0 : if(out.shape()(0) < (x/2+1))
126 0 : throw(AipsError("out shape has to be x/2+1 in size for real to complex FFT2D"));
127 0 : for(uInt k=1; k < shp.nelements(); ++k){
128 0 : if(shp(k) != out.shape()(k))
129 0 : throw(AipsError("shapes of out lattice does not match in lattice for FFT2D"));
130 : }
131 0 : Long numplanes=shp.product()/x/y;
132 0 : IPosition blc(shp.nelements(), 0);
133 0 : IPosition shape=in.shape();
134 :
135 0 : for (uInt ax=2; ax < shp.nelements(); ++ax)
136 0 : shape(ax)=1;
137 0 : IPosition outshape=shape;
138 0 : outshape(0)=x/2+1;
139 :
140 0 : Array<Complex> arr;
141 0 : Array<Float> arrf;
142 : Bool isRef;
143 : Bool del;
144 : Bool delf;
145 : Complex *scr;
146 : Float *scrf;
147 :
148 :
149 :
150 0 : for (Long n=0; n< numplanes; ++n){
151 0 : isRef=out.getSlice(arr, blc, outshape);
152 0 : scr=arr.getStorage(del);
153 : ///Use this method rather than arrf=in.getSlice(blc,shape)
154 : ///as this may be a reference ..the other is a copy always...
155 : /// can gain 0.8s or so for a 10000x10000 array circa 2016
156 0 : in.getSlice( arrf, blc, shape);
157 0 : scrf=arrf.getStorage(delf);
158 0 : r2cFFT(scr, scrf, x, y);
159 0 : arr.putStorage(scr, del);
160 0 : arrf.putStorage(scrf, delf);
161 :
162 0 : if(!isRef){
163 0 : out.putSlice(arr, blc);
164 :
165 : }
166 : //Now calculate next plane
167 :
168 0 : Bool addNextAx=true;
169 0 : for (uInt ax=2; ax < shp.nelements(); ++ax){
170 0 : if(addNextAx){
171 0 : blc(ax) +=1;
172 0 : addNextAx=false;
173 : }
174 0 : if(blc(ax)> shp(ax)-1){
175 0 : blc(ax)=0;
176 0 : addNextAx=true;
177 : }
178 :
179 : }
180 :
181 : }
182 0 : }
183 0 : void FFT2D::r2cFFT(Complex*& out, Float*& in, Long x, Long y){
184 0 : if(x%2 != 0 || y%2 != 0)
185 0 : throw(AipsError("Programmer error: FFT2D does not deal with odd numbers on x-y plane"));
186 0 : fftShift(in, x, y);
187 0 : doFFT(out, in, x, y);
188 : //fft1_p.plan_r2c(IPosition(2,x,y), in, out);
189 : //fft1_p.r2c(IPosition(2,x,y), in, out);
190 : //flipArray out is of shape x/2+1, y
191 0 : Complex* scr=out;
192 0 : Matrix<Complex> tmpo(x/2+1, y/2);
193 : Bool gool;
194 0 : Complex* tmpptr=tmpo.getStorage(gool);
195 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
196 : for (Long jj=0; jj< y/2; ++jj){
197 : for(Long ii=0; ii < (x/2+1); ++ii){
198 : tmpptr[jj*(x/2+1)+ii]=scr[jj*(x/2+1)+ii];
199 : scr[jj*(x/2+1)+ii]=scr[(y/2)*(x/2+1)+jj*(x/2+1)+ii];
200 : }
201 : }
202 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
203 : for (Long jj=0; jj< y/2; ++jj){
204 : for(Long ii=0; ii < x/2; ++ii){
205 : scr[(y/2)*(x/2+1)+jj*(x/2+1)+ii]=tmpptr[jj*(x/2+1)+ii];
206 : }
207 : }
208 :
209 0 : }
210 4 : void FFT2D::c2cFFT(Lattice<Complex>& inout, Bool toFreq){
211 4 : IPosition shp=inout.shape();
212 4 : if(shp.nelements() <2)
213 0 : throw(AipsError("Lattice has to be 2 dimensional to use FFT2D"));
214 4 : Long x= inout.shape()(0);
215 4 : Long y=inout.shape()(1);
216 4 : Long numplanes=inout.shape().product()/x/y;
217 4 : IPosition blc(inout.shape().nelements(), 0);
218 4 : IPosition shape=inout.shape();
219 12 : for (uInt ax=2; ax < shp.nelements(); ++ax)
220 8 : shape(ax)=1;
221 4 : Array<Complex> arr;
222 : Bool isRef;
223 : Bool del;
224 : Complex *scr;
225 :
226 8 : for (Long n=0; n< numplanes; ++n){
227 4 : isRef=inout.getSlice(arr, blc, shape);
228 4 : scr=arr.getStorage(del);
229 4 : c2cFFT(scr, x, y, toFreq);
230 4 : arr.putStorage(scr, del);
231 4 : if(!isRef)
232 0 : inout.putSlice(arr, blc);
233 : //Now calculate next plane
234 4 : Bool addNextAx=true;
235 12 : for (uInt ax=2; ax < shp.nelements(); ++ax){
236 8 : if(addNextAx){
237 8 : blc(ax) +=1;
238 8 : addNextAx=false;
239 : }
240 8 : if(blc(ax)> shp(ax)-1){
241 8 : blc(ax)=0;
242 8 : addNextAx=true;
243 : }
244 :
245 : }
246 : }
247 4 : }
248 0 : void FFT2D::c2cFFTInDouble(Lattice<Complex>& inout, Bool toFreq){
249 0 : IPosition shp=inout.shape();
250 0 : if(shp.nelements() <2)
251 0 : throw(AipsError("Lattice has to be 2 dimensional to use FFT2D"));
252 0 : Long x= inout.shape()(0);
253 0 : Long y=inout.shape()(1);
254 0 : Long numplanes=inout.shape().product()/x/y;
255 0 : IPosition blc(inout.shape().nelements(), 0);
256 0 : IPosition shape=inout.shape();
257 0 : for (uInt ax=2; ax < shp.nelements(); ++ax)
258 0 : shape(ax)=1;
259 0 : Array<Complex> arr;
260 : Bool isRef;
261 : Bool del, delD;
262 : Complex *scr;
263 : DComplex *scrD;
264 0 : Array<DComplex> arrD(shape);
265 0 : scrD=arrD.getStorage(delD);
266 0 : for (Long n=0; n< numplanes; ++n){
267 0 : isRef=inout.getSlice(arr, blc, shape);
268 0 : scr=arr.getStorage(del);
269 0 : complexConvert(scrD, scr, (ooLong)shape(0)*(ooLong)shape(1), False);
270 0 : c2cFFT(scrD, x, y, toFreq);
271 0 : complexConvert(scrD, scr, (ooLong)shape(0)*(ooLong)shape(1), True);
272 0 : arr.putStorage(scr, del);
273 0 : if(!isRef)
274 0 : inout.putSlice(arr, blc);
275 : //Now calculate next plane
276 0 : Bool addNextAx=true;
277 0 : for (uInt ax=2; ax < shp.nelements(); ++ax){
278 0 : if(addNextAx){
279 0 : blc(ax) +=1;
280 0 : addNextAx=false;
281 : }
282 0 : if(blc(ax)> shp(ax)-1){
283 0 : blc(ax)=0;
284 0 : addNextAx=true;
285 : }
286 :
287 : }
288 : }
289 0 : arrD.putStorage(scrD, delD);
290 0 : }
291 0 : void FFT2D::c2cFFT(Lattice<DComplex>& inout, Bool toFreq){
292 0 : IPosition shp=inout.shape();
293 0 : if(shp.nelements() <2)
294 0 : throw(AipsError("Lattice has to be 2 dimensional to use FFT2D"));
295 0 : Long x= inout.shape()(0);
296 0 : Long y=inout.shape()(1);
297 0 : Long numplanes=inout.shape().product()/x/y;
298 0 : IPosition blc(inout.shape().nelements(), 0);
299 0 : IPosition shape=inout.shape();
300 0 : for (uInt ax=2; ax < shp.nelements(); ++ax)
301 0 : shape(ax)=1;
302 0 : Array<DComplex> arr;
303 : Bool isRef;
304 : Bool del;
305 : DComplex *scr;
306 :
307 0 : for (Long n=0; n< numplanes; ++n){
308 0 : isRef=inout.getSlice(arr, blc, shape);
309 0 : scr=arr.getStorage(del);
310 0 : c2cFFT(scr, x, y, toFreq);
311 0 : arr.putStorage(scr, del);
312 0 : if(!isRef)
313 0 : inout.putSlice(arr, blc);
314 : //Now calculate next plane
315 0 : Bool addNextAx=true;
316 0 : for (uInt ax=2; ax < shp.nelements(); ++ax){
317 0 : if(addNextAx){
318 0 : blc(ax) +=1;
319 0 : addNextAx=false;
320 : }
321 0 : if(blc(ax)> shp(ax)-1){
322 0 : blc(ax)=0;
323 0 : addNextAx=true;
324 : }
325 :
326 : }
327 : }
328 0 : }
329 :
330 4 : void FFT2D::c2cFFT(Complex*& out, Long x, Long y, Bool toFreq){
331 4 : if(x%2 != 0 || y%2 !=0)
332 0 : throw(AipsError("Programmer error: FFT2D does not deal with odd numbers on x-y plane"));
333 4 : fftShift(out, x, y, true);
334 4 : doFFT(out, x, y, toFreq);
335 : /*Int dim[2]={Int(x), Int(y)};
336 : if(toFreq){
337 :
338 : planC2C_p=fftwf_plan_dft(2, dim, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out), FFTW_FORWARD, FFTW_ESTIMATE);
339 :
340 : //fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
341 : }
342 : else{
343 : planC2C_p=fftwf_plan_dft(2, dim, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out), FFTW_BACKWARD, FFTW_ESTIMATE);
344 : // fft1_p.plan_c2c_backward(IPosition(2, x, y), out);
345 : }
346 : fftwf_execute(planC2C_p);
347 : */
348 4 : fftShift(out, x, y, toFreq);
349 :
350 4 : }
351 0 : void FFT2D::c2cFFT(DComplex*& out, Long x, Long y, Bool toFreq){
352 0 : if(x%2 != 0 || y%2 !=0)
353 0 : throw(AipsError("Programmer error: FFT2D does not deal with odd numbers on x-y plane"));
354 0 : fftShift(out, x, y, true);
355 0 : doFFT(out, x, y, toFreq);
356 0 : fftShift(out, x, y, toFreq);
357 :
358 0 : }
359 0 : void FFT2D::doFFT(DComplex*& out, Long x, Long y, Bool toFreq){
360 0 : if(useFFTW_p){
361 : //Will need to seperate the plan from the execute if we want to run this in multiple threads
362 0 : Int dim[2]={Int(y), Int(x)};
363 0 : if(toFreq){
364 0 : if(!planC2CD_forw_p){
365 0 : planC2CD_forw_p=fftw_plan_dft(2, dim, reinterpret_cast<fftw_complex *>(out), reinterpret_cast<fftw_complex *>(out), FFTW_FORWARD, FFTW_ESTIMATE);
366 0 : fftw_execute(planC2CD_forw_p);
367 0 : nx_p=x;
368 0 : ny_p=y;
369 : }
370 : else{
371 0 : if((nx_p != x) || (ny_p !=y)) {
372 0 : throw_programmer_error(nx_p, ny_p, x, y, __FILE__, __LINE__);
373 : }
374 0 : fftw_execute_dft(planC2CD_forw_p, reinterpret_cast<fftw_complex *>(out), reinterpret_cast<fftw_complex *>(out));
375 : }
376 : //fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
377 : }
378 : else{
379 0 : if(!planC2CD_back_p){
380 0 : planC2CD_back_p=fftw_plan_dft(2, dim, reinterpret_cast<fftw_complex *>(out), reinterpret_cast<fftw_complex *>(out), FFTW_BACKWARD, FFTW_ESTIMATE);
381 0 : fftw_execute(planC2CD_back_p);
382 0 : nx_p=x;
383 0 : ny_p=y;
384 : }
385 : else{
386 0 : if((nx_p != x) || (ny_p !=y)) {
387 0 : throw_programmer_error(nx_p, ny_p, x, y, __FILE__, __LINE__);
388 : }
389 0 : fftw_execute_dft(planC2CD_back_p, reinterpret_cast<fftw_complex *>(out), reinterpret_cast<fftw_complex *>(out));
390 : }
391 : // fft1_p.plan_c2c_backward(IPosition(2, x, y), out);
392 : }
393 : }
394 : else{
395 0 : throw(AipsError("Double precision FFT with FFTPack is not implemented"));
396 : }
397 0 : }
398 4 : void FFT2D::doFFT(Complex*& out, Long x, Long y, Bool toFreq){
399 4 : if(useFFTW_p){
400 : //Will need to seperate the plan from the execute if we want to run this in multiple threads
401 4 : Int dim[2]={Int(y), Int(x)};
402 4 : if(toFreq){
403 2 : if(!planC2C_forw_p){
404 2 : planC2C_forw_p=fftwf_plan_dft(2, dim, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out), FFTW_FORWARD, FFTW_ESTIMATE);
405 2 : fftwf_execute(planC2C_forw_p);
406 2 : nx_p=x;
407 2 : ny_p=y;
408 :
409 : }
410 : else{
411 0 : if((nx_p != x) || (ny_p !=y)) {
412 0 : throw_programmer_error(nx_p, ny_p, x, y, __FILE__, __LINE__);
413 : }
414 0 : fftwf_execute_dft(planC2C_forw_p, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out) );
415 : }
416 : //fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
417 : }
418 : else{
419 2 : if(!planC2C_back_p){
420 2 : planC2C_back_p=fftwf_plan_dft(2, dim, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out), FFTW_BACKWARD, FFTW_ESTIMATE);
421 2 : fftwf_execute(planC2C_back_p);
422 2 : nx_p=x;
423 2 : ny_p=y;
424 : }
425 : else{
426 0 : if((nx_p != x) || (ny_p !=y)) {
427 0 : throw_programmer_error(nx_p, ny_p, x, y, __FILE__, __LINE__);
428 : }
429 0 : fftwf_execute_dft(planC2C_back_p, reinterpret_cast<fftwf_complex *>(out), reinterpret_cast<fftwf_complex *>(out) );
430 : }
431 : // fft1_p.plan_c2c_backward(IPosition(2, x, y), out);
432 : }
433 :
434 :
435 : }
436 : else{
437 : Int ier;
438 0 : Int x1=Int(x);
439 0 : Int y1=Int(y);
440 0 : if(wsave_p.size()==0){
441 0 : wsave_p.resize(2*x1*y1+15);
442 0 : lsav_p=2*x1*y1+15;
443 0 : Float *wsaveptr=wsave_p.data();
444 0 : FFTPack::cfft2i(x1, y1, wsaveptr, lsav_p, ier);
445 : }
446 0 : std::vector<Float> work(2*x1*y1);
447 0 : Int lenwrk=2*x1*y1;
448 0 : Float* workptr=work.data();
449 0 : Float* wsaveptr=wsave_p.data();
450 0 : if(toFreq)
451 0 : FFTPack::cfft2f(x1, x1, y1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
452 : else
453 0 : FFTPack::cfft2b(x1, x1, y1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
454 0 : }
455 4 : }
456 0 : void FFT2D::doFFT(Complex*& out, Float*& in, Long x, Long y){
457 0 : if(useFFTW_p){
458 0 : Int dim[2]={Int(y), Int(x)};
459 0 : if(!planR2C_p){
460 0 : planR2C_p=fftwf_plan_dft_r2c(2, dim, in, reinterpret_cast<fftwf_complex *>(out), FFTW_ESTIMATE);
461 :
462 : //fft1_p.plan_c2c_forward(IPosition(2, x, y), out);
463 :
464 0 : fftwf_execute(planR2C_p);
465 0 : nx_p=x;
466 0 : ny_p=y;
467 : }
468 : else{
469 0 : if((nx_p != x) || (ny_p !=y)) {
470 0 : throw_programmer_error(nx_p, ny_p, x, y, __FILE__, __LINE__);
471 : }
472 0 : fftwf_execute_dft_r2c(planR2C_p, in, reinterpret_cast<fftwf_complex *>(out));
473 : }
474 : }
475 : else{
476 : /*
477 : Int ier;
478 : Int x1=Int(x);
479 : Int y1=Int(y);
480 : if(wsave_p.size()==0){
481 : wsave_p.resize(2*x1*y1+15);
482 : lsav_p=2*x1*y1+15;
483 : Float *wsaveptr=wsave_p.data();
484 : FFTPack::cfft2i(x1, y1, wsaveptr, lsav_p, ier);
485 : }
486 : std::vector<Float> work(2*x1*y1);
487 : Int lenwrk=2*x1*y1;
488 : Float* workptr=work.data();
489 : Float* wsaveptr=wsave_p.data();
490 : if(toFreq)
491 : FFTPack::cfft2f(x1, y1, x1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
492 : else
493 : FFTPack::cfft2b(x1, y1, x1, out, wsaveptr, lsav_p, workptr, lenwrk, ier);
494 : */
495 0 : throw(AipsError("Not implemented FFTPack r2c yet"));
496 : }
497 0 : }
498 :
499 8 : void FFT2D::fftShift(Complex*& s, Long x, Long y, Bool toFreq){
500 : ////Lets try our own flip
501 :
502 : Bool gool;
503 8 : Complex* scr=s;
504 : {
505 8 : Matrix<Complex> tmpo(x/2, y/2);
506 8 : Complex* tmpptr=tmpo.getStorage(gool);
507 : ////TEST
508 : //omp_set_num_threads(1);
509 : /////
510 : /*
511 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
512 : for (Long jj=0; jj< y/2; ++jj){
513 : for(Long ii=0; ii < x/2; ++ii){
514 : tmpptr[jj*x/2+ii]=scr[(y/2-jj-1)*x+(x/2-ii-1)];
515 : scr[(y/2)*x+(jj*x+x/2)+ii]=scr[jj*x+ii];
516 : }
517 : }
518 : */
519 8 : Float divid=1.0f;
520 8 : if(!toFreq)
521 2 : divid=1.0f/(Float(x)*Float(y));
522 8 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr, divid)
523 : for (Long jj=0; jj< y/2; ++jj){
524 : for(Long ii=0; ii < x/2; ++ii){
525 : tmpptr[jj*x/2+ii]=scr[(y/2)*x+(jj*x+x/2)+ii]*divid;
526 : scr[(y/2)*x+(jj*x+x/2)+ii]=scr[jj*x+ii]*divid;
527 : }
528 : }
529 8 : #pragma omp parallel for default(none) firstprivate(x,y, tmpptr, scr)
530 : for (Long jj=0; jj< y/2; ++jj){
531 : for(Long ii=0; ii < x/2; ++ii){
532 : scr[jj*x+ii]=tmpptr[jj*x/2+ii];
533 : }
534 : }
535 8 : #pragma omp parallel for default(none) firstprivate(x,y, tmpptr, scr, divid)
536 : for (Long jj=0; jj< y/2; ++jj){
537 : for(Long ii=0; ii < x/2; ++ii){
538 : tmpptr[jj*x/2+ii]=scr[(jj*x+x/2)+ii]*divid;
539 : scr[(jj*x+x/2)+ii]=scr[(y/2)*x+jj*x+ii]*divid;
540 : }
541 : }
542 8 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
543 : for (Long jj=0; jj< y/2; ++jj){
544 : for(Long ii=0; ii < x/2; ++ii){
545 : scr[(y/2)*x+jj*x+ii]=tmpptr[jj*x/2+ii];
546 : }
547 : }
548 8 : tmpo.putStorage(tmpptr, gool);
549 8 : }
550 :
551 : ////
552 :
553 : //if(rot)
554 : /*{
555 :
556 : Matrix<Complex> tmpo(x, y/2);
557 : Complex* tmpptr=tmpo.getStorage(gool);
558 : for (Long jj=0; jj< y/2; ++jj){
559 : for(Long ii=0; ii < x; ++ii){
560 : tmpptr[jj*x+ii]=scr[(y-jj-1)*x+(x-ii-1)];
561 : scr[(y-jj-1)*x+(x-ii-1)]=scr[jj*x+ii];
562 : }
563 : }
564 : for (Long jj=0; jj< y/2; ++jj){
565 : for(Long ii=0; ii < x; ++ii){
566 : scr[jj*x+ii]= tmpptr[jj*x+ii];
567 : }
568 : }
569 : }*/
570 :
571 :
572 :
573 :
574 8 : }
575 :
576 0 : void FFT2D::fftShift(DComplex*& s, Long x, Long y, Bool toFreq){
577 : ////Lets try our own flip
578 :
579 : Bool gool;
580 0 : DComplex* scr=s;
581 : {
582 0 : Matrix<DComplex> tmpo(x/2, y/2);
583 0 : DComplex* tmpptr=tmpo.getStorage(gool);
584 : ////TEST
585 : //omp_set_num_threads(1);
586 : /////
587 : /*
588 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
589 : for (Long jj=0; jj< y/2; ++jj){
590 : for(Long ii=0; ii < x/2; ++ii){
591 : tmpptr[jj*x/2+ii]=scr[(y/2-jj-1)*x+(x/2-ii-1)];
592 : scr[(y/2)*x+(jj*x+x/2)+ii]=scr[jj*x+ii];
593 : }
594 : }
595 : */
596 0 : Double divid=1.0f;
597 0 : if(!toFreq)
598 0 : divid=1.0f/(Double(x)*Double(y));
599 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr, divid)
600 : for (Long jj=0; jj< y/2; ++jj){
601 : for(Long ii=0; ii < x/2; ++ii){
602 : tmpptr[jj*x/2+ii]=scr[(y/2)*x+(jj*x+x/2)+ii]*divid;
603 : scr[(y/2)*x+(jj*x+x/2)+ii]=scr[jj*x+ii]*divid;
604 : }
605 : }
606 0 : #pragma omp parallel for default(none) firstprivate(x,y, tmpptr, scr)
607 : for (Long jj=0; jj< y/2; ++jj){
608 : for(Long ii=0; ii < x/2; ++ii){
609 : scr[jj*x+ii]=tmpptr[jj*x/2+ii];
610 : }
611 : }
612 0 : #pragma omp parallel for default(none) firstprivate(x,y, tmpptr, scr, divid)
613 : for (Long jj=0; jj< y/2; ++jj){
614 : for(Long ii=0; ii < x/2; ++ii){
615 : tmpptr[jj*x/2+ii]=scr[(jj*x+x/2)+ii]*divid;
616 : scr[(jj*x+x/2)+ii]=scr[(y/2)*x+jj*x+ii]*divid;
617 : }
618 : }
619 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
620 : for (Long jj=0; jj< y/2; ++jj){
621 : for(Long ii=0; ii < x/2; ++ii){
622 : scr[(y/2)*x+jj*x+ii]=tmpptr[jj*x/2+ii];
623 : }
624 : }
625 0 : tmpo.putStorage(tmpptr, gool);
626 0 : }
627 :
628 : ////
629 :
630 : //if(rot)
631 : /*{
632 :
633 : Matrix<Complex> tmpo(x, y/2);
634 : Complex* tmpptr=tmpo.getStorage(gool);
635 : for (Long jj=0; jj< y/2; ++jj){
636 : for(Long ii=0; ii < x; ++ii){
637 : tmpptr[jj*x+ii]=scr[(y-jj-1)*x+(x-ii-1)];
638 : scr[(y-jj-1)*x+(x-ii-1)]=scr[jj*x+ii];
639 : }
640 : }
641 : for (Long jj=0; jj< y/2; ++jj){
642 : for(Long ii=0; ii < x; ++ii){
643 : scr[jj*x+ii]= tmpptr[jj*x+ii];
644 : }
645 : }
646 : }*/
647 :
648 :
649 :
650 :
651 0 : }
652 0 : void FFT2D::fftShift(Float*& s, Long x, Long y){
653 : ////Lets try our own flip
654 :
655 : Bool gool;
656 0 : Float* scr=s;
657 0 : Matrix<Float> tmpo(x/2, y/2);
658 0 : Float* tmpptr=tmpo.getStorage(gool);
659 : ////TEST
660 : //omp_set_num_threads(1);
661 : /////
662 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
663 : for (Long jj=0; jj< y/2; ++jj){
664 : for(Long ii=0; ii < x/2; ++ii){
665 : tmpptr[jj*x/2+ii]=scr[(y/2)*x+(jj*x+x/2)+ii];
666 : scr[(y/2)*x+(jj*x+x/2)+ii]=scr[jj*x+ii];
667 : }
668 : }
669 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
670 : for (Long jj=0; jj< y/2; ++jj){
671 : for(Long ii=0; ii < x/2; ++ii){
672 : scr[jj*x+ii]=tmpptr[jj*x/2+ii];
673 : }
674 : }
675 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
676 : for (Long jj=0; jj< y/2; ++jj){
677 : for(Long ii=0; ii < x/2; ++ii){
678 : tmpptr[jj*x/2+ii]=scr[(jj*x+x/2)+ii];
679 : scr[(jj*x+x/2)+ii]=scr[(y/2)*x+jj*x+ii];
680 : }
681 : }
682 0 : #pragma omp parallel for default(none) firstprivate(x, y, tmpptr, scr)
683 : for (Long jj=0; jj< y/2; ++jj){
684 : for(Long ii=0; ii < x/2; ++ii){
685 : scr[(y/2)*x+jj*x+ii]=tmpptr[jj*x/2+ii];
686 : }
687 : }
688 :
689 :
690 :
691 0 : }
692 0 : void FFT2D::complexConvert(DComplex*& scrD, Complex*& scr, const ooLong len, const Bool down){
693 0 : if(down){
694 0 : for(ooLong k=0; k< len; ++k){
695 0 : scr[k]=(Complex)(scrD[k]);
696 : }
697 : }
698 : else{
699 0 : for(ooLong k=0; k< len; ++k){
700 0 : scrD[k]=(DComplex)(scr[k]);
701 : }
702 : }
703 0 : }
704 0 : std::tuple<Long, Long> FFT2D::getShape() {
705 :
706 0 : return std::make_tuple(nx_p, ny_p);
707 : }
708 :
709 :
710 : } //# NAMESPACE CASA - END
|
