2 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
6 * A PICTURE RECONSTRUCTION PROGRAM *
8 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
10 snfft.cpp,v 1.2 2008/08/25 16:11:09 jklukowska Exp
14 PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX
15 THREE-DIMENSIONAL ARRAY.
17 CALL SNFFT(X,N,INVDIR)
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
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
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.
42 FFT2P, FTODD, UNKPS, FREE
45 INTEGER FUNCTION MAVIRT(I,J,Y)
48 THIS FUNCTION IS ONLY REQUIRED DURING LOADING,BUT IS NOT
49 CALLED DURING EXECUTION.
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.
57 FOR INVDIR = +1 THE FORIER TRANSFORM OF COMPLEX ARRAY
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
64 WHERE WI IS THE N(I)-TH ROOT OF UNITY AND L1=P*K,
66 FOR INVDIR = -1 THE INVERSE TRANSFORM OF COMPLEX ARRAY
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
75 WHERE NORM = N(1)*N(2)*N(3).
76 ..................................................................
91 void snfft(REAL* x, INTEGER* n, INTEGER invdir)
113 REAL* iftaux = NULL; //wei 3/2005
121 fprintf(output, "\n **** parameter invdir of function snfft() should be +1 or -1");
122 fprintf(output, "\n **** program aborted\n");
126 for(i = 0; i < 3; i++) {
129 fprintf(output, "\n **** parameter n[%d] of function snfft() should be positive", i);
130 fprintf(output, "\n **** program aborted\n");
136 upfac = n[0] * n[1] * n[2];
137 for(i = 0; i < 3; i++) {
143 dostep = nxtelt * ifft;
144 uppbnd = dostep * (upfac - 1) + 1;
146 ipf = new INTEGER[7];
148 primfc(ifft, &m, &nipf, ipf, 7);
150 ipivot = ifft / (1 << (2 * (m / 2)));
154 ipusiz = (3 * ipivot) / 2 + 2;
159 ndftau = 2 * ipf[nipf - 1] * nxtelt;
160 itsize = 2 * ipf[nipf - 1];
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;
174 for(j = 0; j < uppbnd; j += dostep) {
176 fft2p(nxtelt, nxtelt, iftaux, x+ind1, m, -invdir, ix, nipf, ipf, ipu, ivspt, it, FALSE);
178 // test /////////////////////////////////////////
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]);
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]);
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]);
201 /////////////////////////////////////////////////
204 //////////////////////////////
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]);
212 /////////////////////////////////////////////////////
216 delete[] ipf; // bug 92 - Lajos - 03/02/2005
223 lhi = 2 * n[0] * n[1] * n[2];
225 facnor = ((REAL) 2.0) / ((REAL) lhi);
227 for(l = 0; l < lhi; l++) {
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