src/snark/ftodd.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/ftodd.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  THE SUBROUTINES FFT2P , FTODD , UNKPS ARE FAST FOUTIER ROUTINES.
  11  THE SUBROUTINE SNFFT ,DESCRIBED IN THE SNARK05 MANUAL , MAKES
  12  USE OF THESE SUBROUTINES.
  13  THE SUBROUTINES HAVE BEEN TAKEN FROM
  14  POLGE,R.J., AND BHAGAVAN,B.K., EFFICIENT FAST FOURIER TRANSFORM
  15  PROGRAMS FOR ARBITRARY FACTORS WITH ONE STEP LOOP UNSCRAMBLING,
  16  IEEE TRANSACTIONS ON COMPUTERS, MAY 1976.
  17  */
  18
  19 #include <cstdio>
  20 #include <cmath>
  21
  22 #include "blkdta.h"
  23
  24 #include "ftodd.h"
  25
  26 void ftodd(
  27 INTEGER imax,
  28 INTEGER istep,
  29 REAL* xs,
  30 REAL* x,
  31 INTEGER n,
  32 INTEGER* nk,
  33 REAL* v,
  34 REAL tpis,
  35 INTEGER nipf,
  36 INTEGER* ipf,
  37 INTEGER mc,
  38 INTEGER k4m4,
  39 INTEGER nipfx,
  40 INTEGER* ipfx,
  41 INTEGER* ipu,
  42 INTEGER* ix,
  43 INTEGER* nfixp,
  44 REAL* vspt,
  45 REAL* t,
  46 BOOLEAN virt)
  47 {
  48         REAL vs[2];
  49
  50         INTEGER i;
  51         INTEGER ip;
  52         INTEGER sk;
  53         REAL fix;
  54         INTEGER kt;
  55         REAL temp;
  56         INTEGER kss;
  57         INTEGER ks;
  58         INTEGER kp;
  59         INTEGER kgp;
  60         INTEGER ktas;
  61         INTEGER kta;
  62         INTEGER ia;
  63         INTEGER l;
  64         INTEGER ind2;
  65         INTEGER ind;
  66         INTEGER indx;
  67         INTEGER kpm1;
  68         INTEGER jt;
  69         INTEGER ktip;
  70         INTEGER k;
  71         INTEGER jlim;
  72         INTEGER kradx;
  73         INTEGER kradxm;
  74         INTEGER j;
  75         INTEGER ir;
  76         INTEGER inc;
  77         INTEGER jl;
  78         INTEGER irm1;
  79         INTEGER kiu;
  80         INTEGER kix;
  81         INTEGER inp;
  82         INTEGER in;
  83
  84         for (i = 0; i < nipf; i++)
  85         {
  86                 ip = ipf[i];
  87                 sk = *nk / ip;
  88                 CA(t, 0, 0)= 1.0;
  89                 CA(t, 1, 0)= 0.0;
  90
  91                 fix = tpis / ip;
  92
  93                 vs[0] = (REAL) cos(fix);
  94                 vs[1] = (REAL) sin(fix);
  95
  96                 CA(vspt, 0, 0)= 1.0;
  97                 CA(vspt, 1, 0)= 0.;
  98
  99                 for (kt = 1; kt < ip; kt++)
 100                 {
 101                         temp = CA(vspt, 0, kt-1)* vs[0] - CA(vspt, 1, kt-1) * vs[1];
 102                         CA(vspt, 1, kt) = CA(vspt, 0, kt-1) * vs[1] + CA(vspt, 1, kt-1) * vs[0];
 103                         CA(vspt, 0, kt) = temp;
 104                 }
 105
 106                 vs[0] = 1.;
 107                 vs[1] = 0.;
 108
 109                 for (kss = 0; kss < sk; kss++)
 110                 {
 111
 112                         ks = kss;
 113
 114                         for (kp = 1; kp < ip; kp++)
 115                         {
 116                                 temp = CA(t, 0, kp-1)* vs[0] - CA(t, 1, kp-1) * vs[1];
 117                                 CA(t, 1, kp) = CA(t, 0, kp-1) * vs[1] + CA(t, 1, kp-1) * vs[0];
 118                                 CA(t, 0, kp) = temp;
 119                         }
 120
 121                         for (kgp = 0; kgp < n; kgp += *nk)
 122                         {
 123                                 ktas = ks + kgp;
 124                                 kta = ktas;
 125
 126                                 for (kt = 0; kt < ip; kt++)
 127                                 {
 128                                         ia = kta * istep;
 129
 130                                         for (l = 0; l < imax; l++)
 131                                         {
 132                                                 ind2 = l + ia;
 133                                                 ind = l + kt * istep;
 134                                                 CA(xs, 0, ind)= CA(x, 0, ind2);
 135                                                 CA(xs, 1, ind)= CA(x, 1, ind2);
 136                                         }
 137                                         kta += sk;
 138                                 }
 139
 140                                 ia = ktas * istep;
 141
 142                                 for (l = 0; l < imax; l++)
 143                                 {
 144                                         ind = l + ia;
 145                                         CA(x, 0, ind)= CA(xs, 0, l);
 146                                         CA(x, 1, ind)= CA(xs, 1, l);
 147                                 }
 148
 149                                 for (kp = 1; kp < ip; kp++)
 150                                 {
 151                                         for (l = 0; l < imax; l++)
 152                                         {
 153                                                 ind = l + kp * istep;
 154                                                 indx = l + ia;
 155                                                 CA(x, 0, indx)+= CA(xs, 0, ind);
 156                                                 CA(x, 1, indx)+= CA(xs, 1, ind);
 157                                         }
 158                                 }
 159
 160                                 kta = ktas;
 161                                 for (kp = 1; kp < ip; kp++)
 162                                 {
 163                                         kta += sk;
 164                                         ia = kta * istep;
 165
 166                                         for (l = 0; l < imax; l++)
 167                                         {
 168                                                 ind = l + ia;
 169                                                 CA(x, 0, ind)= CA(xs, 0, l);
 170                                                 CA(x, 1, ind)= CA(xs, 1, l);
 171                                         }
 172
 173                                         kpm1 = kp;
 174                                         jt = kpm1;
 175
 176                                         for (kt = 1; kt < ip; kt++)
 177                                         {
 178                                                 jt = jt - ip * (jt / ip); // jt %= ip // ??? mk
 179                                                 for (l = 0; l < imax; l++)
 180                                                 {
 181                                                         ind = l + kt * istep;
 182                                                         indx = l + ia;
 183                                                         CA(x, 0, indx)+= CA(vspt, 0, jt) * CA(xs, 0, ind) - CA(vspt, 1, jt) * CA(xs, 1, ind);
 184                                                         CA(x, 1, indx)+= CA(vspt, 0, jt) * CA(xs, 1, ind) + CA(vspt, 1, jt) * CA(xs, 0, ind);
 185                                                 }
 186                                                 jt += kpm1;
 187                                         }
 188
 189                                         for (l = 0; l < imax; l++)
 190                                         {
 191                                                 indx = l + ia;
 192                                                 temp = CA(x, 0, indx)* CA(t, 0, kp) - CA(x, 1, indx) * CA(t, 1, kp);
 193                                                 CA(x, 1, indx)= CA(x, 0, indx) * CA(t, 1, kp) + CA(x, 1, indx) * CA(t, 0, kp);
 194                                                 CA(x, 0, indx)= temp;
 195                                         }
 196                                 }
 197                         }
 198                         temp = v[0] * vs[0] - v[1] * vs[1];
 199                         vs[1] = v[0] * vs[1] + v[1] * vs[0];
 200                         vs[0] = temp;
 201                 }
 202                 *nk = sk;
 203
 204                 CA(xs, 0, 0)= v[0];
 205                 CA(xs, 1, 0)= v[1];
 206
 207                 for (l = 1; l < ip; l++)
 208                 {
 209                         temp = CA(xs, 0, 0)* v[0] - CA(xs, 1, 0) * v[1];
 210                         CA(xs, 1, 0) = CA(xs, 0, 0) * v[1] + CA(xs, 1, 0) * v[0];
 211                         CA(xs, 0, 0) = temp;
 212                 }
 213
 214                 v[0] = CA(xs, 0, 0);
 215                 v[1] = CA(xs, 1, 0);
 216         }
 217
 218         if (nipfx < 2)
 219                 return;
 220
 221         ktip = mc - mc / 2;
 222
 223         ipfx[0] = 4;
 224         if (mc == 1)
 225                 ipfx[0] = 2;
 226         if (ktip == 2)
 227                 ipfx[1] = 2;
 228
 229         for (k = 0; k < nipf; k++)
 230         {
 231                 ind = k + ktip;
 232                 ipfx[ind] = ipf[k];
 233         }
 234
 235         ix[0] = 0;
 236         jlim = 1;
 237
 238         for (i = 0; i < nipfx; i++)
 239         {
 240                 kradx = ipfx[i];
 241                 kradxm = kradx - 1;
 242
 243                 for (j = 0; j < jlim; j++)
 244                 {
 245                         ir = j;
 246                         inc = 0;
 247                         ix[j] *= kradx;
 248
 249                         for (jl = 0; jl < kradxm; jl++)
 250                         {
 251                                 ir += jlim;
 252                                 inc++;
 253                                 ix[ir] = ix[j] + inc;
 254                         }
 255                 }
 256                 jlim *= kradx;
 257         }
 258
 259         irm1 = ir;
 260         kiu = 0;
 261
 262         for (kix = 1; kix < irm1; kix++)
 263         {
 264                 if (!(ix[kix] != kix))
 265                 {
 266                         (*nfixp)++;
 267                         ipu[kiu] = ix[kix] * k4m4;
 268                         ix[kix] = -ix[kix];
 269                         kiu++;
 270                 }
 271         }
 272
 273         for (kix = 1; kix < irm1; kix++)
 274         {
 275                 if (!(ix[kix] < 0))
 276                 {
 277                         inp = kix;
 278                         do
 279                         {
 280                                 in = ix[inp];
 281                                 ix[inp] = -ix[inp];
 282                                 inp = in;
 283                                 ipu[kiu] = in * k4m4;
 284
 285                                 kiu++;
 286
 287                         } while (inp != kix);
 288
 289                         ipu[kiu] = 0;
 290                         kiu++;
 291                 }
 292         }
 293
 294         ipu[kiu] = n;
 295         return;
 296 }