src/snark/fft2p.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/fft2p.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  fft2p.cpp,v 1.2 2008/08/25 16:11:08 jklukowska Exp
  11
  12  THE SUBROUTINES FFT2P , FTODD , UNKPS ARE FAST FOUTIER ROUTINES.
  13  THE SUBROUTINE SNFFT ,DESCRIBED IN THE SNARK05 MANUAL , MAKES
  14  USE OF THESE SUBROUTINES.
  15  THE SUBROUTINES HAVE BEEN TAKEN FROM
  16  POLGE,R.J., AND BHAGAVAN,B.K., EFFICIENT FAST FOURIER TRANSFORM
  17  PROGRAMS FOR ARBITRARY FACTORS WITH ONE STEP LOOP UNSCRAMBLING,
  18  IEEE TRANSACTIONS ON COMPUTERS, MAY 1976.
  19  */
  20
  21 #include <cstdio>
  22 #include <cmath>
  23
  24 #include "blkdta.h"
  25 #include "ftodd.h"
  26 #include "unkps.h"
  27 #include "consts.h"
  28 #include "fft2p.h"
  29
  30 void fft2p(
  31 INTEGER imax,
  32 INTEGER istep,
  33 REAL* xs,
  34 REAL* x,
  35 INTEGER m,
  36 INTEGER invdir,
  37 INTEGER* ix,
  38 INTEGER nipf,
  39 INTEGER* ipf,
  40 INTEGER* ipu,
  41 REAL* vspt,
  42 REAL* t,
  43 BOOLEAN virt)
  44 {
  45
  46         REAL v[2], vg[2];
  47
  48         //BOOLEAN dovirt;
  49         INTEGER sk, ipfx[8];
  50
  51         REAL tpis;
  52         INTEGER m2;
  53         INTEGER mc;
  54         INTEGER irp;
  55         INTEGER k;
  56         INTEGER nfixp;
  57         INTEGER nipfx;
  58         INTEGER ir;
  59         INTEGER n;
  60         INTEGER k4m4;
  61         INTEGER m1;
  62         INTEGER mx;
  63         REAL fi;
  64         INTEGER nk;
  65         INTEGER j;
  66         INTEGER kp;
  67         INTEGER ia;
  68         INTEGER ib;
  69         INTEGER l;
  70         INTEGER index1;
  71         INTEGER index2;
  72         REAL temp;
  73
  74         //dovirt = virt && (istep > 1);
  75
  76         tpis = ((REAL) 2.0) * Consts.pi / ((REAL) invdir);
  77         m2 = m / 2;
  78         mc = m - 2 * m2;
  79         irp = 1;
  80
  81         if (nipf != 0)
  82         {
  83                 for (k = 0; k < nipf; k++)
  84                 {
  85                         irp *= ipf[k];
  86                 }
  87                 nfixp = 0;
  88         }
  89
  90         nipfx = nipf + mc;
  91         ir = irp * (1 << mc);
  92         n = (1 << (m - mc)) * ir;
  93
  94         k4m4 = 1 << m2;
  95
  96         m1 = m2 + mc;
  97         mx = m;
  98         fi = tpis / n;
  99
 100         v[0] = (REAL) cos(fi);
 101         v[1] = (REAL) sin(fi);
 102
 103         nk = n;
 104         if (mx != 0)
 105         {
 106
 107                 L1: sk = nk / 2;
 108                 vg[0] = 1.0;
 109                 vg[1] = 0.0;
 110                 for (j = 0; j < sk; j++)
 111                 {
 112                         for (k = j; k < n; k += nk)
 113                         {
 114                                 kp = k + sk;
 115                                 ia = k * istep;
 116
 117                                 ib = kp * istep;
 118
 119                                 for (l = 0; l < imax; l++)
 120                                 {
 121
 122                                         index1 = l + ia;
 123                                         index2 = l + ib;
 124
 125                                         CA(xs, 0, 0)= CA(x, 0, index1);
 126                                         CA(xs, 1, 0)= CA(x, 1, index1);
 127
 128                                         CA(x, 0, index1)+= CA(x, 0, index2);
 129                                         CA(x, 1, index1)+= CA(x, 1, index2);
 130
 131                                         CA(xs, 0, 0)-= CA(x, 0, index2);
 132                                         CA(xs, 1, 0)-= CA(x, 1, index2);
 133
 134                                         CA(x, 0, index2)= vg[0] * CA(xs, 0, 0) - vg[1] * CA(xs, 1, 0);
 135                                         CA(x, 1, index2)= vg[0] * CA(xs, 1, 0) + vg[1] * CA(xs, 0, 0);
 136                                 }
 137                         }
 138                         temp = vg[0] * v[0] - vg[1] * v[1];
 139                         vg[1] = vg[0] * v[1] + vg[1] * v[0];
 140                         vg[0] = temp;
 141                 }
 142                 nk = sk;
 143                 temp = (v[0] - v[1]) * (v[0] + v[1]);
 144                 v[1] = ((REAL) 2.0) * v[0] * v[1];
 145                 v[0] = temp;
 146                 m1--;
 147                 mx--;
 148                 if (m1 > 0)
 149                         goto L1
 150                         ;
 151                 if (mx != 0)
 152                 {
 153                         m1 = m2;
 154                 }
 155         }
 156
 157         if (irp != 1)
 158         {
 159
 160                 ftodd(imax, istep, xs, x, n, &nk, v, tpis, nipf, ipf, mc, k4m4, nipfx, ipfx, ipu, ix, &nfixp, vspt, t, virt);
 161
 162                 irp = 1;
 163         }
 164
 165         if (m1 != 0)
 166                 goto L1;
 167
 168         unkps(imax, istep, xs, x, m, nipfx, 2, mc, ir, n, ix, ipu, nfixp, k4m4, virt);
 169         return;
 170 }