src/snark/snfft.cpp

   1 /*
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                           *
   4  *                              S N A R K   1 4                              *
   5  *                                                                           *
   6  *                     A PICTURE RECONSTRUCTION PROGRAM                      *
   7  *                                                                           *
   8  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   9
  10  snfft.cpp,v 1.2 2008/08/25 16:11:09 jklukowska Exp
  11
  12     SUBROUTINE SNFFT
  13     PURPOSE
  14        PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX
  15        THREE-DIMENSIONAL ARRAY.
  16     USAGE
  17        CALL SNFFT(X,N,INVDIR)
  18     DESCRIPTION
  19        X     - AS INPUT, X CONTAINS THE COMPLEX THREE-DIMENSIONAL
  20                ARRAY TO BE TRANSFORMED. THE REAL PART OF X(K,L,M)
  21                IS STORED IN VECTOR FASHION IN A CELL WITH INDEX
  22                2*(K*N(3)*N(2)+L*N(3)+M)+1, THE IMAGINARY PART IS
  23                STORED IN THE CELL IMMEDIATELY FOLLOWING. NOTE THAT
  24                THE SUBSCRIPT K INCREASES MOST RAPIDLY, AND M LEAST
  25                RAPIDLY.
  26                AS OUTPUT, X CONTAINS THE COMPLEX FOURIER TRANSFORM
  27                OR THE INVERSE TRANSFORM, DEPENDING ON INVDIR.
  28        N     - A THREE CELL VECTOR WHICH DETERMINES THE SIZE OF
  29                THE THREE DIMENSIONS OF X.
  30        INVDIR- INTEGER, WHICH SHOULD BE EITHER +1 OR -1, WHICH
  31                DETERMINES WHETHER THE DIRECT OR INVERSE TRANSFORM
  32                IS TO BE PERFORMED.
  33     REMARKS
  34        THIS SUBROUTINE CAN BE USED IN A SNARK-ENVIRONMENT ONLY]
  35        THE WORKSPACE THIS SUBROUTINE REQUIRES IS ALLOCATED FROM
  36        THE BLANK COMMON BLOCK, USING THE SNARK SUBROUTINE ALLOC.
  37        THIS WORKSPACE IS RELEASED UPON EXIT OF SNFFT.
  38        THE SUBROUTINE CAN BE USED FOR ANY VALUES OF THE VECTOR N.
  39        IT IS POSSIBLE THAT THE SUBROUTINE CAUSES A HALT , DUE TO
  40        THE FACT THAT NOT ENOUGH SPACE IN BLANK COMMON IS AVAILABLE.
  41     SUBROUTINES REQUIRED
  42        FFT2P, FTODD, UNKPS, FREE
  43     FUNCTIONS REQUIRED
  44        ALLOC, MAVIRT
  45        INTEGER FUNCTION MAVIRT(I,J,Y)
  46        INTEGER I,J
  47        REAL Y(1)
  48        THIS FUNCTION IS ONLY REQUIRED DURING LOADING,BUT IS NOT
  49        CALLED DURING EXECUTION.
  50     SOURCE
  51        THE SUBROUTINES FFT2P, FTODD AND UNKPS ARE FROM
  52           POLGE,R.J. AND BHAGAVAN,B.K.,EFFICIENT FAST FOURIER
  53                TRANSFORM PROGRAMS FOR ARBITRARY FACTORS WITH
  54                ONE STEP LOOP UNSCRAMBLING. IEEE TRANSACTIONS ON
  55                COMPUTERS, MAY 1976, PP. 534-539.
  56     METHOD
  57        FOR INVDIR = +1 THE FORIER TRANSFORM OF COMPLEX ARRAY
  58        X IS OBTAINED.
  59
  60            N(1)-1  N(2)-1  N(3)-1            -L1  -L2  -L3
  61   X(K,L,M)=SUM     SUM     SUM     X(P,Q,R)*W1  *W2  *W3
  62            P=0     Q=0     R=0
  63
  64        WHERE WI IS THE N(I)-TH ROOT OF UNITY AND L1=P*K,
  65        L2=L*Q AND L3=R*M.
  66        FOR INVDIR = -1 THE INVERSE TRANSFORM  OF COMPLEX ARRAY
  67        X IS OBTAINED.
  68
  69  X(K,L,M) =
  70
  71     1      N(1)-1  N(2)-1  N(3)-1             L1   L2   L3
  72    ---- *  SUM     SUM     SUM     X(P,Q,R)*W1  *W2  *W3
  73    NORM    P=0     Q=0     R=0
  74
  75        WHERE NORM = N(1)*N(2)*N(3).
  76  ..................................................................
  77 */
  78
  79 #include <cstdlib>
  80 #include <cstdio>
  81 #include <cmath>
  82
  83 #include "blkdta.h"
  84 #include "uiod.h"
  85
  86 #include "primfc.h"
  87 #include "fft2p.h"
  88 #include "snfft.h"
  89
  90
  91 void snfft(REAL* x, INTEGER* n, INTEGER invdir)
  92 {
  93   INTEGER i;
  94   INTEGER nxtelt;
  95   INTEGER upfac;
  96   INTEGER ifft;
  97   INTEGER dostep;
  98   INTEGER uppbnd;
  99   INTEGER* ipf;
 100   INTEGER m;
 101   INTEGER nipf;
 102   INTEGER ipivot;
 103   INTEGER ndftau;
 104   INTEGER ipusiz;
 105   INTEGER itsize;
 106   //INTEGER* ind;
 107   INTEGER ind1;
 108   INTEGER ixsiz;
 109   REAL* ivspt = NULL;
 110   REAL* it = NULL;
 111   INTEGER* ix = NULL;
 112   INTEGER* ipu = NULL;
 113   REAL* iftaux = NULL;    //wei 3/2005
 114   INTEGER j;
 115   INTEGER lhi;
 116   REAL facnor;
 117   INTEGER l;
 118
 119   if(abs(invdir) != 1)
 120   {
 121     fprintf(output, "\n **** parameter invdir of function snfft() should be +1 or -1");
 122     fprintf(output, "\n **** program aborted\n");
 123     exit(0);
 124   }
 125
 126   for(i = 0; i < 3; i++) {
 127     if(n[i] < 1)
 128     {
 129       fprintf(output, "\n **** parameter n[%d] of function snfft() should be positive", i);
 130       fprintf(output, "\n **** program aborted\n");
 131       exit(0);
 132     }
 133   }
 134
 135   nxtelt = 1;
 136   upfac = n[0] * n[1] * n[2];
 137   for(i = 0; i < 3; i++) {
 138
 139     ifft = n[2-i];
 140
 141     if(ifft != 1) {
 142       upfac /= ifft;
 143       dostep = nxtelt * ifft;
 144       uppbnd = dostep * (upfac - 1) + 1;
 145
 146       ipf = new INTEGER[7];
 147
 148       primfc(ifft, &m, &nipf, ipf, 7);
 149
 150       ipivot = ifft / (1 << (2 * (m / 2)));
 151
 152       ndftau = 2 * nxtelt;
 153
 154       ipusiz = (3 * ipivot) / 2 + 2;
 155
 156       itsize = 2;
 157
 158       if(nipf != 0) {
 159         ndftau = 2 * ipf[nipf - 1] * nxtelt;
 160         itsize = 2 * ipf[nipf - 1];
 161       }
 162       ixsiz = MAX0(ipusiz, (1 << (m/2)));
 163       ivspt = new REAL[itsize];
 164       it = new REAL[itsize];
 165       ix = new INTEGER[ixsiz];
 166       ipu = new INTEGER[ipusiz];
 167       iftaux = new REAL[ndftau];
 168       for (j=0;j<itsize;j++) ivspt[j]=0;
 169       for (j=0;j<itsize;j++) it[j]=0;
 170       for (j=0;j<ixsiz;j++) ix[j]=0;
 171       for (j=0;j<ipusiz;j++) ipu[j]=0;
 172       for (j=0;j<ndftau;j++) iftaux[j]=0;
 173
 174       for(j = 0; j < uppbnd; j += dostep) {
 175         ind1 = 2 * j;
 176         fft2p(nxtelt, nxtelt, iftaux, x+ind1, m, -invdir, ix, nipf, ipf, ipu, ivspt, it, FALSE);
 177
 178         // test /////////////////////////////////////////
 179         /*
 180         fprintf(output,"\n snfft ivspt");
 181         for(int xx = 0; xx < itsize; xx++) {
 182           if((xx % 7) == 0) fprintf(output,"\n");
 183           fprintf(output," %9.4f", ivspt[xx]);
 184         }
 185         */
 186         /*
 187         fprintf(output,"\n snfft it");
 188         for(int xx = 0; xx < itsize; xx++) {
 189           if((xx % 7) == 0) fprintf(output,"\n");
 190           fprintf(output," %9.4f", it[xx]);
 191         }
 192         */
 193
 194         /*
 195         fprintf(output,"\n snfft x");
 196         for(int xx = 0; xx < 2 * n[0] * n[1] * n[2]; xx++) {
 197           if((xx % 7) == 0) fprintf(output,"\n");
 198           fprintf(output," %9.4f", x[xx]);
 199         }
 200         */
 201         /////////////////////////////////////////////////
 202       }
 203
 204       //////////////////////////////
 205       /*
 206         fprintf(output,"\n snfft x");
 207         for(int xx = 0; xx < 2 * n[0] * n[1] * n[2]; xx++) {
 208           if((xx % 7) == 0) fprintf(output,"\n");
 209           fprintf(output," %9.4f", x[xx]);
 210         }
 211       */
 212       /////////////////////////////////////////////////////
 213
 214       nxtelt *= ifft;
 215
 216       delete[] ipf;  // bug 92 - Lajos - 03/02/2005
 217
 218     }
 219   }
 220
 221   if(invdir == -1) {
 222
 223     lhi = 2 * n[0] * n[1] * n[2];
 224
 225     facnor = ((REAL) 2.0) / ((REAL) lhi);
 226
 227     for(l = 0; l < lhi; l++) {
 228       x[l] *= facnor;
 229     }
 230   }
 231
 232
 233        if(ivspt != NULL)  delete[] ivspt;
 234        if(it != NULL)  delete[] it;
 235        if(ix != NULL)  delete[] ix;
 236        if(ipu != NULL)  delete[] ipu;
 237        if(iftaux != NULL)  delete[] iftaux;       //wei 3/2005
 238 }