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/rtfort.cpp $
5 $LastChangedRevision: 95 $
6 $Date: 2014-07-02 19:43:36 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 ************************************************************
12 SUBROUTINE RTFORT(X, N, INVDIR)
15 RTFORT COMPUTES THE HERMITIAN FOURIER SPECTRUM OF A ONE OR TWO
16 DIMENSIONAL ARRAY OF REAL DATA (IF INVDIR=1), OR THE REAL DATA
17 CORRESPONDING TO THIS SPECTRUM (IF INVDIR=-1).
20 TO COMPUTE DIRECT FOURIER TRANSFORM OF REAL ARRAY X (INVDIR=1)
21 CALL RTFORT( X, N, INVDIR )
22 HALF OF THE CONJUGATE-SYMMETRIC (HERMITIAN) SPECTRUM
23 IS NOW STORED IN X AS COMPLEX FOURIER COEFFICIENTS.
24 TO COMPUTE INVERSE TRANSFORM OF SEMI-SPECTRUM (INVDIR = -1)
25 CALL RTFORT(X, N, INVDIR)
26 ARRAY X NOW CONTAINS REAL DATA.
30 ONE DIMENSIONAL, N(1)=1 : X CONTAINS 2*N(2) REAL NUMBERS IN
31 CONSECUTIVE LOCATIONS AS INPUT TO DIRECT TRANSFORM OR AS RESULT
32 OF INVERSE TRANSFORM. THE SEMI-SPECTRUM CONSISTS OF ( N(2) + 1
33 COMPLEX COEFFICIENTS STORED IN 2*N(2) + 2 LOCATIONS OF X.
35 TWO DIMENSIONAL : X IS A REAL ARRAY OF (FAST*SLOW) DIMENSIONS
36 2N(2) * N(1) AS INPUT TO DIRECT TRANSFORM OR AS RESULT OF
37 INVERSE TRANSFORM. THE SEMI-SPECTRUM CONSISTS OF
38 ( N(2) + 1 ) * N(1) (FAST*SLOW) COMPLEX FOURIER COEFFICIENTS
39 STORED IN ( 2*N(2) + 2 ) * N(1) LOCATIONS OF X.
41 N.B. X MUST BE DIMENSIONED TO ACCOMODATE THE SEMI-SPECTRUM,
42 WHICH REQUIRES MORE STORAGE THAN ITS CORRESPONDING REAL ARRAY.
45 A 3-ELEMENT VECTOR, N(1)= SLOW DIMENSION OF X ARRAY
46 2*N(2) = FAST DIMENSION OF REAL X ARRAY (AN EVEN NUMBER)
47 N(3) = 1 (ONLY 2-DIMENSIONAL ARRAYS ALLOWED)
48 THE VALUES OF N(1) AND N(2) MUST BE COMPATIBLE WITH THE
49 PARTICULAR COMPLEX FFT ROUTINE USED (E.G. SNFFT, HARM).
51 INVDIR = 1 FOR DIRECT (FORWARD) TRANSFORM
52 = -1 FOR INVERSE TRANSFORM
54 SUBROUTINE CALLED: TRANSM
56 BASED ON R2FORT SUBROUTINE WRITTEN BY T.M.PETERS, JUNE 1971,
57 FOR USE WITH RADIX TWO FFT. RTFORT EXTENDS CAPABILITY TO
58 RADIX ODD TRANSFORMS. R.M. LEWITT, SEPT 1977.
60 **************************************************************
73 void rtfort(REAL* x, INTEGER* nin, INTEGER invdir)
119 // INITIALISE PARAMETERS FOR RECURSIVE GENERATION OF SINES
121 sd = Consts.pid2 / m;
122 cd = ((REAL) 2.0) * SQR((REAL) sin(sd));
124 sd = (REAL) -sin(sd + sd);
126 // ASSUMING THE DIRECT FOURIER TRANSFORM HAS EXP(-I...)
127 // IF DIRECT FT ROUTINE HAS EXP(+I...) FACTOR, THEN
128 // REPLACE ABOVE STATEMENT BY
135 // ONE DIMENSIONAL CASE
142 // REPEAT FIRST TWO REAL NUMBERS AT ARRAY END FOR PERIODICITY.
154 // APPLY TRANSFORM SYMMETRY RELATIONS
155 for (if_1 = 0; if_1 < m1; if_1 += 2)
158 transm(x, if_1, ib, sn, cn);
159 ct = cn - (cd * cn + sd * sn);
160 sn = sn + (sd * cn - cd * sn);
161 ca = ((REAL) 0.5) / (SQR(ct) + SQR(sn)) + ((REAL) 0.5);
169 // FOR INVERSE TRANSFORM, ZERO LAST TWO ARRAY ELEMENTS
176 // TWO DIMENSIONAL CASE
177 // MTYPE = 0 (1) FOR M ODD (EVEN)
178 mtype = (m / 2) * 2 - m + 1;
183 // FORWARD TRANSFORM. EXPAND ARRAY TO TRANSFORM FORMAT
185 for (k = 1; k <= ns; k++)
187 ishift = (ns - k) * 2;
188 for (i = 0; i < mm; i++)
191 x[index + ishift] = x[index];
193 igap = index + ishift + mm;
195 x[igap + 1] = x[index + 1];
207 for (i1 = 0; i1 < m; i1 += 2)
209 // K=0 CASE SEPARATELY
212 transm(x, j1, j2, sn, cn);
215 for (k1 = 1; k1 < ns; k1++)
220 transm(x, j1, j2, sn, cn);
223 ct = cn - (cd * cn + sd * sn);
224 sn = sn + (sd * cn - cd * sn);
225 ca = ((REAL) 0.5) / (SQR(ct) + SQR(sn)) + ((REAL) 0.5);
233 // M EVEN, TAKE I1=M1 AS SPECIAL CASE
234 m1 -= 1; // decrement m1 - swr 11/25/04
236 transm(x, m1, m1, sn, cn);
239 for (k1 = 2; k1 <= nl; k1++)
240 { // k => k1 - swr 11/25/04
243 j2 = jt - j1 - 2; // decrement by 2 - swr 11/25/04
244 transm(x, j1, j2, sn, cn);
251 // INVERSE TRANSFORM, REVERT TO ORIGINAL FORMAT
253 for (k = 1; k < ns; k++)
256 for (i = 1; i <= mm; i++)
259 x[index] = x[index + ishift];
264 // ZERO REMAINDER OF REAL ARRAY
267 for (i = j1; i < j2; i++)