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) $
8 ***********************************************************
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.
84 for (i = 0; i < nipf; i++)
93 vs[0] = (REAL) cos(fix);
94 vs[1] = (REAL) sin(fix);
99 for (kt = 1; kt < ip; kt++)
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;
109 for (kss = 0; kss < sk; kss++)
114 for (kp = 1; kp < ip; kp++)
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];
121 for (kgp = 0; kgp < n; kgp += *nk)
126 for (kt = 0; kt < ip; kt++)
130 for (l = 0; l < imax; l++)
133 ind = l + kt * istep;
134 CA(xs, 0, ind)= CA(x, 0, ind2);
135 CA(xs, 1, ind)= CA(x, 1, ind2);
142 for (l = 0; l < imax; l++)
145 CA(x, 0, ind)= CA(xs, 0, l);
146 CA(x, 1, ind)= CA(xs, 1, l);
149 for (kp = 1; kp < ip; kp++)
151 for (l = 0; l < imax; l++)
153 ind = l + kp * istep;
155 CA(x, 0, indx)+= CA(xs, 0, ind);
156 CA(x, 1, indx)+= CA(xs, 1, ind);
161 for (kp = 1; kp < ip; kp++)
166 for (l = 0; l < imax; l++)
169 CA(x, 0, ind)= CA(xs, 0, l);
170 CA(x, 1, ind)= CA(xs, 1, l);
176 for (kt = 1; kt < ip; kt++)
178 jt = jt - ip * (jt / ip); // jt %= ip // ??? mk
179 for (l = 0; l < imax; l++)
181 ind = l + kt * istep;
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);
189 for (l = 0; l < imax; l++)
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;
198 temp = v[0] * vs[0] - v[1] * vs[1];
199 vs[1] = v[0] * vs[1] + v[1] * vs[0];
207 for (l = 1; l < ip; l++)
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];
229 for (k = 0; k < nipf; k++)
238 for (i = 0; i < nipfx; i++)
243 for (j = 0; j < jlim; j++)
249 for (jl = 0; jl < kradxm; jl++)
253 ix[ir] = ix[j] + inc;
262 for (kix = 1; kix < irm1; kix++)
264 if (!(ix[kix] != kix))
267 ipu[kiu] = ix[kix] * k4m4;
273 for (kix = 1; kix < irm1; kix++)
283 ipu[kiu] = in * k4m4;
287 } while (inp != kix);