File: fft_utils.i

package info (click to toggle)
yorick-yutils 1.3.0-2
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 604 kB
  • ctags: 17
  • sloc: makefile: 105; python: 12
file content (876 lines) | stat: -rw-r--r-- 32,256 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
/*
 * fft_utils.i --
 *
 *	Useful routines for FFT operations in Yorick.
 *
 * Copyright (c) 1995-2003 Eric THIEBAUT.
 *
 * This program is free software; you can redistribute it and/or  modify it
 * under the terms of the GNU General Public License  as  published  by the
 * Free Software Foundation; either version 2 of the License,  or  (at your
 * option) any later version.
 *
 * This program is distributed in the hope  that  it  will  be  useful, but
 * WITHOUT  ANY   WARRANTY;   without   even   the   implied   warranty  of
 * MERCHANTABILITY or  FITNESS  FOR  A  PARTICULAR  PURPOSE.   See  the GNU
 * General Public License for more details (to receive a  copy  of  the GNU
 * General Public License, write to the Free Software Foundation, Inc., 675
 * Mass Ave, Cambridge, MA 02139, USA).
 *
 * History:
 *	$Id: fft_utils.i,v 1.1.1.1 2007/12/11 23:55:14 frigaut Exp $
 *	$Log: fft_utils.i,v $
 *	Revision 1.1.1.1  2007/12/11 23:55:14  frigaut
 *	Initial Import - yorick-yutils
 *	
 *	Revision 1.9  2004/10/11 11:18:21  eric
 *	 -  Fix fft_dist() so that it is as flexible as, e.g., array() for
 *	    the dimension list.
 *
 *	Revision 1.8  2004/08/31 16:20:22  eric
 *	 - New routines for Fourier interpolation: fft_fine_shift,
 *	   fft_unphasor, fft_interp, fft_interp_complex, fft_interp_real.
 *
 *	Revision 1.7  2003/08/23 09:55:57  eric
 *	 - new functions: fft_gaussian_mtf and fft_gaussian_psf;
 *	 - *** POSSIBLE INCOMPATIBILITY *** change order of args in
 *	   fft_recenter and use fft_setup to speed up FFT's;
 *
 *	Revision 1.6  2003/01/31 15:58:32  eric
 *	 - Added new routines: fft_recenter, fft_smooth and reverse_dims.
 *
 *	Revision 1.5  2002/11/20 09:20:56  eric
 *	 - new keywords SQUARE and NYQUIST in fft_dist
 *	 - make use of grow_dimlist in "utils.i"
 *
 *	Revision 1.4  2002/11/14 10:59:38  eric
 *	 - new graphics routines: fft_plh, fft_plc, fft_plfc
 *	 - new routines: abs2, fft_roll_1d, fft_roll_2d, __fft & Co
 *	 - introduction documentation by: help, fft_utils;
 *
 *	Revision 1.3  2001/04/23 14:16:15  eric
 *	New routine: fft_symmetric_index.
 *
 *	Revision 1.2  2001/04/06 16:42:59  eric
 *	 - new routine: fft_plg
 *	 - fixed routine: fft_recenter_at_max
 *
 *	Revision 1.1  2001/03/23 16:45:54  eric
 *	Initial revision
 *
 *-----------------------------------------------------------------------------
 */

require, "utils.i";

local fft_utils;
/* DOCUMENT: FFT utility routines in "fft_utils.i"
     This package  is mainly written  to deal with the  particular indexing
     rules in FFT transformed arrays.  The following routines are provided:
       abs2         - squared absolute value.
       fft_best_dim - get best dimension to compute FFT.
       fft_centroid - get centroid in FFT arrays.
       fft_convolve - compute discrete convolution thanks to FFT.
       fft_dist     - compute length of FFT frequencies/coordinates.
       fft_fine_shift, fft_unphasor - shift/roll an array by non-integer
                      offset by means of Fourier interpolation.
       fft_gaussian_mtf - compute Gaussian modulation transfer function.
       fft_gaussian_psf - compute Gaussian point spread function.
       fft_indgen   - generate index of FFT frequencies/coordinates.
       fft_interp, fft_interp_complex, fft_interp_real - interpolate array
                      at non-integer offsets.
       fft_plc      - plot contours of 2D FFT.
       fft_plfc     - plot filled contours of 2D FFT.
       fft_plg      - plot 1D FFT as curve.
       fft_plh      - plot 1D FFT with stairs.
       fft_pli      - plot 2D FFT as image.
       fft_recenter - recenter array with respect to a template.
       fft_recenter_at_max - recenter FFT arrays at their maximum.
       fft_roll_1d  - roll dimension of 1D arrays.
       fft_roll_2d  - roll dimension of 2D arrays.
       fft_shift_phasor - get complex phasor for arbitrary shift.
       fft_smooth   - smooth an array by convolution with a gaussian.
       fft_symmetric_index - get hermitian-symmetry index for FFT arrays.
       reverse_dims - reverse all dimensions of an array.
       __fft        - expert driver for repeated FFT's with same dimensions.
       __fft_init   - initialization for __fft.

   SEE ALSO: fft, fftw. */

func abs2(x)
/* DOCUMENT abs2(x)
     Returns abs(X)^2
   SEE ALSO: abs. */
{
  if (structof(x) != complex) return x*x;
  y = x.im;
  x = double(x);
  return x*x + y*y;
}

func fft_best_dim(len)
/* DOCUMENT fft_best_dim(len);
     Return the smallest integer which is greater or equal LEN and which is
     a multiple of powers of 2, 3 and/or 5
   SEE ALSO fft_indgen, fft. */
{
  best= 2*len;
  for (i5=1; i5<=len; i5*=5) {
    for (i3= i5; i3<=len; i3*=3) {
      i2= i3;
      while (i2 < len) i2*=2;
      if (i2 == len) return len;
      if (i2-len < best-len) best= i2;
    }
  }
  return best;
}

func fft_indgen(dim) { return (u= indgen(0:dim-1)) - dim*(u > dim/2); }
/* DOCUMENT fft_indgen(len)
     Return FFT frequencies along a dimension of length LEN.
   SEE ALSO: indgen, span, fft_dist, fft_freqlist, fft_symmetric_index. */

func fft_dist(.., nyquist=, square=)
/* DOCUMENT fft_dist(dimlist);
       -or- fft_dist(dim1, dim2, ...);
     Returns Euclidian lenght of spatial frequencies in frequel units for a
     FFT of dimensions DIMLIST.

     If keyword  NYQUIST is true,  the frequel coordinates get  rescaled so
     that the Nyquist frequency is  equal to NYQUIST along every dimension.
     This is obtained by using coordinates:
        (2.0*NYQUIST/DIM(i))*fft_indgen(DIM(i))
     along i-th dimension of lenght DIM(i).

     If  keyword  SQUARE is  true,  the square  of  the  Euclidian norm  is
     returned instead.
     
   SEE ALSO: fft_indgen, fft_symmetric_index. */
{
  /* Build dimension list. */
  local arg, dims;
  while (more_args()) {
    eq_nocopy, arg, next_arg();
    if (is_array(arg) && ((s = structof(arg)) == long || s == int ||
                          s == short || s == char)) {
      /* got an integer array */
      if (! (n = dimsof(arg)(1))) {
        /* got a scalar */
        grow, dims, arg;
      } else if (n == 1 && (n = numberof(arg) - 1) == arg(1)) {
        /* got a vector which is a valid dimension list */
        if (n) grow, dims, arg(2:);
      } else {
        error, "bad dimension list";
      }
    } else if (! is_void(arg)) {
      error, "unexpected data type in dimension list";
    }
  }
  if (is_void(dims)) return 0.0; /* scalar array */
  if (min(dims) <= 0) error, "negative value in dimension list";

  /* Build square radius array one dimension at a time, starting with the
     last dimension. */
  if ((n = numberof(dims)) == 1) return 0.0; /* scalar array */
  if (is_void(nyquist)) {
    r2 = (u = fft_indgen(dims(n)))*u;
    while (--n >= 1) {
      r2 = r2(-,..) + (u = fft_indgen(dims(n)))*u;
    }
  } else {
    nyquist *= 2.0;
    for (k=n ; k>1 ; --k) {
      dim = dims(k);
      u = (nyquist/dim)*fft_indgen(dim);
      r2 = (k<n ? r2(-,..) + u*u : u*u);
    }
  }
  return (square ? r2 : sqrt(r2));
}

func fft_freqlist(dimlist)
/* DOCUMENT ptr = fft_freqlist(dimlist)
     returns a vector of DIMLIST(1) pointers with normalized FFT
     frequencies along all dimensions of DIMLIST (must be
     dimsof(SOME_ARRAY)) with "adequate" geometry:

       *ptr(1) =   (2*pi/dimlist(2))*fft_indgen(dimlist(2));
       *ptr(2) =  [(2*pi/dimlist(3))*fft_indgen(dimlist(3))];
       *ptr(3) = [[(2*pi/dimlist(4))*fft_indgen(dimlist(4))]];
        ...
     
   SEE ALSO: fft_indgen, fft_shift_phasor. */
{
  /* Precompute (scaled) Fourier frequencies along every dimensions (the
     trick is to build these as "vectors" with adequate dimension list so
     as to minimize the number of operations during the search). */
  PI = 3.1415926535897932384626433832795029;
  ndims = dimlist(1);
  ptr = array(pointer, ndims);
  for (k=1 ; k<=ndims ; ++k) {
    len = dimlist(k + 1);
    ws = array(1, k + 1); /* to build dimension list of k-th dimension */
    ws(1) = k;
    ws(0) = len;
    (ws = array(double, ws))(*) = (2.0*PI/len)*fft_indgen(len);
    ptr(k) = &ws;
  }
  return ptr;
}

func fft_smooth(a, fwhm, setup=)
/* DOCUMENT fft_smooth(a, fwhm)
       -or- fft_smooth(a, fwhm, setup=workspace)
     Returns array A smoothed along  all its dimensions by convolution by a
     gaussian  with  full  width  at  half maximum  equals  to  FWHM.   See
     fft_setup for the meaning of keyword SETUP.

   SEE ALSO: fft, fft_setup, fft_gaussian_mtf. */
{
  dims = dimsof(a);
  if (is_void(setup)) setup=fft_setup(dims);
  as_double = (structof(a) != complex);
  a = fft((1.0/numberof(a))*fft_gaussian_mtf(dims, fwhm)*
          fft(a, -1, setup=setup), +1, setup=setup);
  return (as_double ? double(a) : a);
}

func fft_gaussian_psf(dimlist, fwhm) {
  p = exp(-(2.772588722239781237668928485832706272302/fwhm^2)*
          fft_dist(dimlist, square=1)); n = dimsof(p)(1);
  return ((0.9394372786996513337723403284101825868414/fwhm)^n)*p; }
func fft_gaussian_mtf(dimlist, fwhm) {
  nyquist = 0.9433593338763967991867107983208536017193*fwhm;
  return exp(-fft_dist(dimlist, square=1, nyquist=nyquist)); }
/* DOCUMENT fft_gaussian_psf(dimlist, fwhm)
       -or- fft_gaussian_mtf(dimlist, fwhm)
     Returns   normalized   Gaussian  point   spread   function  (PSF)   or
     corresponding modulation  transfer function (MTF)  with dimension list
     DIMLIST  and full  width at  half maximum  equals to  FWHM  along each
     dimensions (in  the PSF space).  Up  to errors due  to limited support
     and/or numerical precision, the PSF and the MTF obey:
     
        sum(PSF) = MTF(1) = 1                       (normalization)
        MTF = fft(PSF, +1)
        PSF = fft(MTF, -1)/numberof(MTF)
        
     where MTF(1) is the 0-th frequency in the MTF.  The standard deviation
     SIGMA and the FWHM are related by:

        FWHM = sqrt(8*log(2))*SIGMA
             ~ 2.354820045031*SIGMA
        
     Note that, owing  to the limited size of  the support and/or numerical
     precision, these properties may not be perfectly met; for that reason,
     _always_ compute directly  what you need, e.g. do not  take the FFT of
     the PSF if what  you need is the MTF.  Also note  that the geometry is
     that of  the FFT and that  for unequal dimension lengths,  the PSF has
     the same width (in "pixels") along every dimension but not the MTF.

   SEE ALSO: fft_dist, fft_smooth. */
/*
 * Notes:
 *   The FFT of a Gaussian of given FWHM is:
 *      exp(-pi^2*fwhm^2*(k/dim)^2/log(16))
 *   where K is the FFT index; hence:
 *      Nyquist = sqrt(pi^2*fwhm^2*(dim/2/dim)^2/4/log(2))
 *              = pi*fwhm/4/sqrt(log(2))
 *              ~ 0.9433593338763967992*fwhm
 *   The Gaussian is (N is the number of dimensions):
 *      ((sqrt(log(16)/pi)/fwhm)^n)*exp(-log(16)*(k/fwhm)^2)
 *     ~ ((0.9394372786996513337/fwhm)^n)*exp(...)
 *
 *   Constants (with 40 significant digits):
 *      pi/4/sqrt(log(2)) = 0.9433593338763967991867107983208536017193;
 *      sqrt(log(16)/pi)  = 0.9394372786996513337723403284101825868414;
 *      log(16)           = 2.772588722239781237668928485832706272302;
 */

/*---------------------------------------------------------------------------*/
/* SIMPLE FAST FOURIER TRANSFORM */

func __fft_init(dimlist)
/* DOCUMENT __fft_init, dimlist;
     Initializes  FFT  workspace  for  further  calls to  __fft  (to  see).
     DIMLIST is the dimension list of the arrays to transform.  The routine
     defines 2 external symbols:     
       __fft_setup  - used to store the FFT workspace)
       __fft_number - used to keep track of the number of calls to __fft
     In order to avoid namespace pollution/clash, a routine that uses __fft
     should declare these symbols as local before calling __fft_init, e.g.:
       local __fft_setup, __fft_number;
       __fft_init, dimlist;

   SEE ALSO: fft_setup, __fft. */
{
  extern __fft_setup, __fft_number;
  dims = dimlist(2:);
  ndims = numberof(dims);
  __fft_setup = array(pointer, ndims);
  __fft_number = 0;
  for (i=1 ; i<=ndims ; ++i) {
    if (! __fft_setup(i)) {
      dim = dims(i);
      ws = array(double, 6*dim + 15);
      fft_init, dim, ws;
      __fft_setup(where(dims == dim)) = &ws;
    }
  }
}

func __fft(x, dir)
/* DOCUMENT __fft(x);
       -or- __fft(x, dir);
     Replacement for Yorick's fft to speed up fast Fourier transforms (FFT)
     in,  e.g.,  iteratives  algorithms  that  necessitate  computation  of
     several  FFT's  of  arrays  with  same dimension  list.   The  FFT  is
     performed on  all dimensions of  X and DIR  must be a  scalar (default
     +1). FFT workspace must be initialized by __fft_init.

   SEE ALSO: __fft_init, fft. */
{
  extern __fft_setup, __fft_number;

  /* Make a private copy of input array even if it is already complex and
     get its dimension list. */
  if (is_void(dir)) dir = +1;
  x = complex(x);
  dims = dimsof(x);
  ndims = dims(1);
  dims = dims(2:0);
  if (structof(__fft_setup) != pointer ||
      numberof(__fft_setup) != ndims) error, "unitialized FFT workspace";

  /* Do the transform along every dimension of X. */
  len = 6*dims + 15; // expected length of workspace vectors
  std = 1;
  top = numberof(x);
  for (i=1 ; i<=ndims ; i++) {
    dim = dims(i);
    ws = __fft_setup(i);
    if (numberof(*ws) != len(i) || structof(*ws) != double)
      error, swrite(format="bad FFT workspace for dimension length %d", dim);
    top /= dim;
    fft_raw, dir, x, std, dim, top, ws;
    std *= dim;
  }

  /* increment counter and return result */
  if (is_void(__fft_number)) __fft_number= 0;
  ++__fft_number;
  return x;
}

/*---------------------------------------------------------------------------*/

func fft_symmetric_index(..)
/* DOCUMENT fft_symmetric_index(dimlist)
       -or- fft_symmetric_index(dim1, dim2, ...);  
     Returns  indices  of  hermitian-symmetry  transform  for  a  FFT  with
     dimension list DIMLIST.  For instance,  if A is a N-dimensional array,
     then:
        AP= A(fft_symmetric_index(dimsof(A)))   
     is  equal to array  A with  its coordinates  negated according  to FFT
     convention:
        AP(X1, X2, ..., XN) = A(-X1, -X2, ..., -XN)
     consequently if A is hermitian then:
        AP= conj(A).
     
   SEE ALSO: fft_indgen. */
{
  /* Build dimension list. */
  dimlist = [0];
  while (more_args()) grow_dimlist, dimlist, next_arg();

  /* Compute result starting by last dimension. */
  local u;
  if ((n = numberof(dimlist)) == 1) return 1; /* scalar array */
  for (k=n ; k>1 ; --k) {
    dim = dimlist(k);
    (q = indgen(dim:1:-1))(1) = 0; // neg. of freq along that dim
    u= k<n ? dim*u(-,..) + q : q;
  }
  return u + 1; /* <== indices start at 1 in Yorick */
}

/*---------------------------------------------------------------------------*/

func _fft_centroid(a1, repeat)
/* DOCUMENT _fft_centroid(a1)
       -or- _fft_centroid(a1, repeat)
     Working routine for fft_centroid: return position of centroid of 1D array
     A1 assuming FFT geometry for the coordinates.

   SEE ALSO fft_centroid. */
{
  dim= numberof(a1);
  u= fft_indgen(dim);
  x0= u(a1(mxx));
  x= double(u);
  do {
    w= x0==0 ? a1 : roll(a1, -x0);
    if (dim%2 == 0) w(dim/2+1)= 0.0; // remove Nyquist frequency
    x1= x0+sum(w*x)/sum(w);
    x0p= x0;
    x0= long(floor(x1+0.5));
  } while (--repeat>0 && abs(x0-x0p)>=1);
  return x1;
}

func fft_centroid(a, repeat)
/* DOCUMENT fft_centroid(a)
       -or- fft_centroid(a, repeat)
     Return  the  position  of  centroid  of N-dimensional  array  A  assuming
     coordinates  along  dimensions  of  A  are  wrapped  as  in  a  FFT  (see
     fft_indgen).  The  algorithm proceeds by computing the  center of gravity
     of A around its  central element which is the maximum of  A for the first
     iteration  and  the  closest  to  the previously  computed  centroid  for
     subsequent  iterations.   The  maximum  number  of  iteration  is  REPEAT
     (default: 3;  in any  cases, at least  one iteration is  performed).  The
     Nyquist frequency along each even dimension is omitted to avoid a bias.

   SEE ALSO fft_indgen. */
{
  if (is_void(repeat)) repeat= 3;
  dims= dimsof(a);
  if ((ndims= dims(1)) <= 2) {
    if (ndims==2) return [_fft_centroid(a(,sum), repeat),
			 _fft_centroid(a(sum,), repeat)];
    if (ndims==1) return _fft_centroid(a, repeat);
    return 0.0; // ndims==0
  }
  result= array(double, ndims);
  for (i=1 ; i<=ndims ; ++i) {
    dim= dims(i+1);
    a1= array(double, numberof(a)/dim, dim);
    a1(*)= (i==ndims ? a(*) : transpose(a, [ndims,i])(*));
    a1= a1(sum,);
    result(i)= _fft_centroid(a1, repeat);
  }
  return result;
}

/*---------------------------------------------------------------------------*/
/* CENTERING / ROLLING */

func reverse_dims(a)
/* DOCUMENT reverse_dims(a)
     Returns array A with all its dimensions reversed.

   SEE ALSO: fft_recenter. */
{
  n = dimsof(a)(1);
  r = ::-1;
  if (n ==  1) return a(r);
  if (n ==  2) return a(r, r);
  if (n ==  3) return a(r, r, r);
  if (n ==  4) return a(r, r, r, r);
  if (n ==  5) return a(r, r, r, r, r);
  if (n ==  6) return a(r, r, r, r, r, r);
  if (n ==  7) return a(r, r, r, r, r, r, r);
  if (n ==  8) return a(r, r, r, r, r, r, r, r);
  if (n ==  9) return a(r, r, r, r, r, r, r, r, r);
  if (n == 10) return a(r, r, r, r, r, r, r, r, r, r);
  if (n == 11) return a(r, r, r, r, r, r, r, r, r, r, r);
  if (n == 12) return a(r, r, r, r, r, r, r, r, r, r, r, r);
  if (n == 13) return a(r, r, r, r, r, r, r, r, r, r, r, r, r);
  if (n == 14) return a(r, r, r, r, r, r, r, r, r, r, r, r, r, r);
  if (n == 15) return a(r, r, r, r, r, r, r, r, r, r, r, r, r, r, r);
  if (n == 16) return a(r, r, r, r, r, r, r, r, r, r, r, r, r, r, r, r);
  if (n == 17) return a(r, r, r, r, r, r, r, r, r, r, r, r, r, r, r, r, r);
  if (n == 18) return a(r, r, r, r, r, r, r, r, r, r, r, r, r, r, r, r, r, r);
  error, "too many dimensions";
}

func fft_recenter(x, template, reverse)
/* DOCUMENT fft_recenter(x, template)
       -or- fft_recenter(x, template, reverse)
     Returns array X rolled so that it matches most closely array TEMPLATE.
     X and TEMPLATE may be arrays  of real or complex numbers but must have
     the same dimension lists.  The  returned value is roll(X, S) where the
     offsets S minimize:

         sum(abs(roll(X, S) - TEMPLATE)^2)

     If optional  argument REVERSE is true,  X is also allowed  to have all
     its dimensions reversed in order to math TEMPLATE.
     
   SEE ALSO: fft, roll, reverse_dims. */
{
  ws = fft_setup(dimsof(a));
  conj_fft_template = conj(fft(template, +1, setup=ws));
  c1 = double(fft(conj_fft_template*fft(x, +1, setup=ws), -1, setup=ws));
  max_c1 = max(c1);
  if (reverse) {
    x2 = reverse_dims(x);
    c2 = double(fft(conj_fft_template*fft(x2, +1, setup=ws), -1, setup=ws));
    if ((max_c2 = max(c2)) > max_c1) {
      eq_nocopy, x, x2;
      i = where(c2 == max_c2);
    } else {
      i = where(c1 == max_c1);
    }
    c2 = [];
  } else {
    i = where(c1 == max_c1);
  }
  c1 = [];
  if (numberof(i) != 1) error, swrite("too many maxima (%d)", numberof(i));
  index = i(1) - 1; /* 0-based position of the maximum of correlation */
  dims = dimsof(x);
  ndims = dims(1);
  dims = dims(2:);
  offset = array(long, ndims);
  for (i=1 ; i<=ndims ; ++i) {
    dim = dims(i);
    offset(i) = index % dim;
    index /= dim;
  }
  //return roll(x, dims - offset);
  return (anyof(offset) ? roll(x, dims - offset) : x);
}

func fft_recenter_at_max(z, middle=)
/* DOCUMENT fft_recenter_at_max(z)
   
     Return Z  rolled so  that its element  with maximum value  (or maximum
     absolute value if  Z is complex) is at the  origin.  If keyword MIDDLE
     is true  (non-zero and non-nil) the  center is at the  middle of every
     dimension  otherwise the  center is  the first  element of  the output
     array (as assumed by the FFT).

   SEE ALSO: roll. */
{
  index = (structof(z) == complex ? abs(z) : z)(*)(mxx) - 1;
  dims = dimsof(z);
  ndims = dims(1);
  dims = dims(2:);
  offset = array(long, ndims);
  for (i=1 ; i<=ndims ; ++i) {
    dim = dims(i);
    offset(i) = index % dim;
    index /= dim;
  }
  if (middle) dims += (dims+1)/2;
  return roll(z, dims - offset);
}

func fft_roll_1d(a, off) {
  if ((dimlist = dimsof(a))(1) != 1) error, "expecting 1D array";
  n = dimlist(2);  k = (n + (off%n))%n; /* wrap offset in the range [0, n-1] */
  if (! k) return a;
  b = array(structof(a), dimlist);
  b(k+1:n) = a(1:n-k); b(1:k) = a(n-k+1:n); return b; }
func fft_roll_2d(a, off1, off2)
/* DOCUMENT fft_roll_1d(v, off)
       -or- fft_roll_2d(m, off1, off2)
     "rolls" dimensions of the vector V (1D array) or matrix M (2D array)
     and return a result with same data type than original array.
     
   SEE ALSO: roll. */
{
  if ((dimlist = dimsof(a))(1) != 2) error, "expecting 2D array";
  n1 = dimlist(2);
  k1 = (n1 + (off1%n1))%n1; /* wrap offset in the range [0, n1-1] */
  n2 = dimlist(3);
  k2 = (n2 + (off2%n2))%n2; /* wrap offset in the range [0, n2-1] */
  case = (! k1) + 2*(! k2);
  if (case == 3) return a;
  b = array(structof(a), dimlist);
  if (case == 0) {
    b((r1=k1+1:n1), (r2=k2+1:n2)) = a((s1=1:n1-k1), (s2=1:n2-k2));
    b(r1, (t2=1:k2)) = a(s1, (u2=n2-k2+1:n2));
    b((t1=1:k1), t2) = a((u1=n1-k1+1:n1), u2);
    b(t1, r2) = a(u1, s2);
  } else if (case == 1) {
    b( , k2+1:n2) = a( , 1:n2-k2);
    b( , 1:k2)    = a( , n2-k2+1:n2);
  } else {
    b(k1+1:n1, ) = a(1:n1-k1, );
    b(1:k1, )    = a(n1-k1+1:n1, );
  }
  return b;
}

/*---------------------------------------------------------------------------*/
/* FOURIER INTERPOLATION */

func fft_shift_phasor(off, u)
/* DOCUMENT fft_shift_phasor(off, dimlist)
       -or- fft_shift_phasor(off, fft_freqlist(dimlist))
     returns complex phasor to apply in FFT space for a shift by OFF cells
     in the real space.  DIMLIST is a list of dimensions -- the second
     callinf sequence is to allow for computing the normalized FFT
     frequencies only once.  The offset OFF must have as many elements as
     PTR or as many as dimensions in DIMLIST (i.e. a shift for each
     dimension) and may be fractionnal.

     This function is intended for Fourier interpolation.  For instance,
     assuming A is 2-D real array:

         z = fft(a, +1);
         u = fft_freqlist(dimsof(a));
         q = fft_unphasor([0.33, -0.47], u);

     Then:

         fft_interp(z, q, real=1);

     yields the value of A interpolated at coordinate (0.33, -0.47) in FFT
     frame, i.e. center of lower left cell is at (0,0).  The shfited version
     of A by (0.33, -0.47) can be obtained by:
     
         fft_shift(z, q, real=1);


   SEE ALSO: fft_freqlist, fft_interp, fft_fine_shift, roll. */
{
  if (structof(u) != pointer) u = fft_freqlist(u);
  ndims = numberof(u);
  for (i=1 ; i<=ndims ; ++i) {
    a = (*u(i))*off(i);
    if (i == 1) p  =  cos(a) + 1i*sin(a);
    else        p *= (cos(a) + 1i*sin(a));
  }
  return p;
}

func fft_unphasor(z, phasor, setup=, real=) {
  z = fft(z*conj(phasor), -1, setup=setup);
  return (1.0/numberof(z))*(real ? double(z) : z); }
func fft_fine_shift(a, off, setup=)
/* DOCUMENT fft_fine_shift(a, off, setup=)
       -or- fft_unphasor(z, phasor, setup=, real=)
     The function fft_fine_shift returns array A shifted by offset OFF
     which can be fractionnal.  Alternatively, the function fft_unphasor
     can be used when the forward FFT of A and/or the complex phasor
     corresponding to the shift are already computed:

       fft_unphasor(fft(A, +1), fft_shift_phasor(OFF, dimsof(A)),
                    real=(structof(A) != complex))

     yields the same result as fft_fine_shift(A, OFF).

     These functions can make use of pre-computed FFT workspace specified
     by keyword SETUP (see fft_setup).
 
     TO-DO: Improve code by not Fourier transforming along direction
            where OFF is zero (or equal to an integer times the length
            of the dimension).
     
   SEE ALSO: fft, fft_setup, fft_shift_phasor, roll. */
{
  real = (structof(a) != complex);
  dims = dimsof(a);
  if (is_void(setup)) setup = fft_setup(dims);
  a = fft(fft(a, +1, setup=setup)*fft_shift_phasor(-off, dims), -1);
  return (1.0/numberof(a))*(real ? double(a) : a);
}

func fft_interp_real(z, phasor)
{ return (1.0/numberof(z))*double(sum(z*phasor)); }
func fft_interp_complex(z, phasor)
{ return (1.0/numberof(z))*sum(z*phasor); }
func fft_interp(a, off, setup=)
/* DOCUMENT fft_interp(a, off, setup=)
       -or- fft_interp_real(z, phasor)
       -or- fft_interp_complex(z, phasor)
     returns value obtained by Fourier interpolation of A at offset OFF.
     The function fft_interp computes the forward FFT of A and can make use
     of pre-computed FFT workspace specified by keyword SETUP.  The two
     other functions (fft_interp_complex, if A is complex; fft_interp_real
     otherwise) are usefull when the forward FFT of A and/or the complex
     phasor corresponding to the shift are already computed, their
     arguments are:

       Z = fft(A, +1);
       PHASOR = fft_shift_phasor(OFF, dimsof(A));
     
   SEE ALSO: fft, fft_fine_shift, fft_shift_phasor. */
{
  real = (structof(a) != complex);
  dims = dimsof(a);
  if (is_void(setup)) setup = fft_setup(dims);
  n = numberof(a);
  a = sum(fft(a, +1, setup=setup)*fft_shift_phasor(-off, dims));
  return (1.0/n)*(real ? double(a) : a);
}

/*---------------------------------------------------------------------------*/
/* GRAPHICS */

func fft_plh(y, scale=, legend=, hide=, type=, width=, color=, smooth=,
             marks=, marker=, mspace=, mphase=)
{ fft_plg, y, scale=scale, legend=legend, hide=hide, type=type, width=width,
    color=color, smooth=smooth, marks=marks, marker=marker, mspace=mspace,
    mphase=mphase, stair=1; }
func fft_plg(y, scale=, legend=, hide=, type=, width=, color=, smooth=,
             marks=, marker=, mspace=, mphase=, stair=)
/* DOCUMENT fft_plg, y;
        -or fft_plh, y;
     Plot 1-D FFT array Y as a curve, taking care of "rolling" Y and setting
     correct coordinates.  Keyword SCALE can be used to indicate the
     "frequel" scale along X-axis (SCALE is a scalar); by default,
     SCALE=1.0.

   KEYWORDS legend, hide, type, width, color, closed, smooth
            marks, marker, mspace, mphase.

   SEE ALSO plh, plg, roll. */
{
  if (is_void(scale)) scale= 1.0;
  else if (! is_array(scale) || dimsof(scale)(1)!=0)
    error, "expecting a scalar for SCALE";
  if (! is_array(y) || (dims= dimsof(y))(1)!=1) error, "expecting 1-D array";

  dim1= dims(2);
  min1= (max1= dim1/2) - dim1 + 1;

  if (stair) {
    // just=
    plh, roll(y, -min1), scale*indgen(min1:max1), legend=legend, hide=hide,
      type=type, width=width, color=color,
      marks=marks, marker=marker, mspace=mspace, mphase=mphase;
  } else {
    plg, roll(y, -min1), scale*indgen(min1:max1), legend=legend, hide=hide,
      type=type, width=width, color=color, smooth=smooth,
      marks=marks, marker=marker, mspace=mspace, mphase=mphase;
  }  
}

func fft_pli(a, scale=, legend=, hide=, top=, cmin=, cmax=)
/* DOCUMENT fft_pli, a;
     Plot 2-D FFT array A as an image, taking care of "rolling" A and setting
     correct world boundaries.  Keyword SCALE can be used to indicate the
     "frequel" scale along both axis (SCALE is a scalar) or along each axis
     (SCALE is a 2-element vector: SCALE=[XSCALE,YSCALE]); by default,
     SCALE=[1.0, 1.0].

   KEYWORDS legend, hide, top, cmin, cmax.

   SEE ALSO pli, fft_roll_2d. */
{
  local scale1, dim1, min1, max1, scale2, dim2, min2, max2;
  __fft_pl2d_limits, a, scale;
  pli, fft_roll_2d(bytscl(a, top=top, cmin=cmin, cmax=cmax), -min1, -min2),
    scale1*(min1 - 0.5), scale2*(min2 - 0.5),
    scale1*(max1 + 0.5), scale2*(max2 + 0.5),
    legend=legend, hide=hide;
}

func fft_plc(a, scale=, levs=, type=, width=, color=, smooth=,
             legend=, hide=, marks=, marker=, mspace=, mphase=)
/* DOCUMENT fft_plc, a;
     Plot contour levels  of a 2-D FFT array A, taking  care of "rolling" A
     and setting  correct world boundaries.   Keyword SCALE can be  used to
     indicate the  "frequel" scale along both  axis (SCALE is  a scalar) or
     along each axis (SCALE  is a 2-element vector: SCALE=[XSCALE,YSCALE]);
     by default, SCALE=[1.0, 1.0].  Other  keywords have same meaning as in
     plc routine.

   KEYWORDS scale, levs, type, width, color, smooth,
            legend, hide, marks, marker, mspace, mphase.

   SEE ALSO plc, roll, fft_plfc. */
{
  local scale1, dim1, min1, max1, scale2, dim2, min2, max2;
  __fft_pl2d_limits, a, scale;
  u1 = (scale1*indgen(min1:max1))(,-:1:dim2);
  u2 = (scale2*indgen(min2:max2))(-:1:dim1,);
  plc, roll(a, [-min1, -min2]), u2, u1,
    levs=levs, type=type, width=width, color=color, smooth=smooth,
    legend=legend, hide=hide,
    marks=marks, marker=marker, mspace=mspace, mphase=mphase;
}

func fft_plfc(a, scale=, levs=, colors=)
/* DOCUMENT fft_plfc, a;
     Plot  filled contour  levels of  a  2-D FFT  array A,  taking care  of
     "rolling" A  and setting correct world boundaries.   Keyword SCALE can
     be used  to indicate the "frequel"  scale along both axis  (SCALE is a
     scalar)   or  along   each  axis   (SCALE  is   a   2-element  vector:
     SCALE=[XSCALE,YSCALE]); by default,  SCALE=[1.0, 1.0].  Other keywords
     have  same meaning  as in  plfc routine.   As with  plfc  routine, the
     actual level values get saved in external symbol plfc_levs.

   KEYWORDS scale, levs, colors.

   SEE ALSO plc, roll, fft_plc. */
{
  local scale1, dim1, min1, max1, scale2, dim2, min2, max2;
  __fft_pl2d_limits, a, scale;
  u1 = (scale1*indgen(min1:max1))(,-:1:dim2);
  u2 = (scale2*indgen(min2:max2))(-:1:dim1,);
  plfc, roll(a, [-min1, -min2]), u2, u1, levs=levs, colors=colors;
}

func __fft_pl2d_limits(z, scale)
/* DOCUMENT __fft_pl2d_limits, z, scale;
     Private routine used by fft_pli, fft_plc and fft_plfc.

   SEE ALSO fft_pli, fft_plc, fft_plfc. */
{
  extern scale1, dim1, min1, max1, scale2, dim2, min2, max2;
  if (is_void(scale)) {
    scale1 = scale2 = 1.0;
  } else if (dimsof(scale)(1) == 0) {
    scale1 = scale2 = scale;
  } else if (numberof(scale) == 2) {
    scale1 = scale(1);
    scale2 = scale(2);
  } else {
    error, "bad number of elements in SCALE";
  }
  if ((dims = dimsof(z))(1) != 2) error, "expecting 2-D array";
  dim1 = dims(2);
  min1 = (max1 = dim1/2) - dim1 + 1;
  dim2 = dims(3);
  min2 = (max2 = dim2/2) - dim2 + 1;
}

/*---------------------------------------------------------------------------*/

func fft_convolve(orig, psf, do_not_roll)
/* DOCUMENT fft_convolve(orig, psf);
       -or- fft_convolve(orig, psf, do_not_roll);
     Return discrete convolution (computed by FFT) of array ORIG by point
     spread function PSF.  Unless argument DO_NOT_ROLL is true, PSF is
     rolled before.  Note: ORIG and PSF must have same dimension list.
     
   SEE ALSO: fft, fft_setup, roll. */
{
  real = (structof(orig) != complex && structof(psf) != complex);
  dims = dimsof(orig);
  ws = fft_setup(dims);
  p = fft(orig, -1, setup=ws) * fft((do_not_roll?psf:roll(psf)), -1, setup=ws);
  orig = psf = []; // possibly free some memory
  fft_inplace, p, +1, setup=ws;
  if (real) p = double(p);
  return (1.0/numberof(p))*p;
}

#if 0
func fft_of_two_real_arrays(a, b, &ft_a, &ft_b, ljdir, rjdir, setup=)
/* DOCUMENT fft_of_two_real_arrays, a, b, ft_a, ft_b, direction;
       -or- fft_of_two_real_arrays, a, b, ft_a, ft_b, ljdir, rjdir;
     Computes the FFT  of arrays A and  B and stores them in  TF_A and FT_B
     respectively.  A and B must have same dimension list.  A single FFT is
     needed.  Agrguments  DIRECTION, LJDIR,  RJDIR, and keyword  SETUP have
     the same meaning as for the fft function (which see).
     
   SEE ALSO: fft_setup, fft_inplace. */
{
  if (structof(a) == complex || structof(b) == complex)
    error, "A and B must be non-complex";
  c = a + 1i*b;
  fft_inplace, c, ljdir, rjdir, setup=setup;
  b = c(fft_symmetric_index(dimsof(c)));
  a = c + b;
  b = c - b;
  ft_a = 0.5*double(a) + 0.5i*(b.im);
  ft_b = 0.5*(a.im)    - 0.5i*double(b);
}
#endif
/*---------------------------------------------------------------------------*/