/* * Copyright(c)1995,97 Mark Olesen * Queen's Univ at Kingston (Canada) * * Permission to use, copy, modify, and distribute this software for * any purpose without fee is hereby granted, provided that this * entire notice is included in all copies of any software which is * or includes a copy or modification of this software and in all * copies of the supporting documentation for such software. * * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR * IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S * UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY * KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS * FITNESS FOR ANY PARTICULAR PURPOSE. * * All of which is to say that you can do what you like with this * source code provided you don't try to sell it as your own and you * include an unaltered copy of this message (including the * copyright). * * It is also implicitly understood that bug fixes and improvements * should make their way back to the general Internet community so * that everyone benefits. * * Changes: * Trivial type modifications by the WebRTC authors. */ /* * File: * WebRtcIsac_Fftn.c * * Public: * WebRtcIsac_Fftn / fftnf (); * * Private: * WebRtcIsac_Fftradix / fftradixf (); * * Descript: * multivariate complex Fourier transform, computed in place * using mixed-radix Fast Fourier Transform algorithm. * * Fortran code by: * RC Singleton, Stanford Research Institute, Sept. 1968 * * translated by f2c (version 19950721). * * int WebRtcIsac_Fftn (int ndim, const int dims[], REAL Re[], REAL Im[], * int iSign, double scaling); * * NDIM = the total number dimensions * DIMS = a vector of array sizes * if NDIM is zero then DIMS must be zero-terminated * * RE and IM hold the real and imaginary components of the data, and return * the resulting real and imaginary Fourier coefficients. Multidimensional * data *must* be allocated contiguously. There is no limit on the number * of dimensions. * * ISIGN = the sign of the complex exponential (ie, forward or inverse FFT) * the magnitude of ISIGN (normally 1) is used to determine the * correct indexing increment (see below). * * SCALING = normalizing constant by which the final result is *divided* * if SCALING == -1, normalize by total dimension of the transform * if SCALING < -1, normalize by the square-root of the total dimension * * example: * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] * * int dims[3] = {n1,n2,n3} * WebRtcIsac_Fftn (3, dims, Re, Im, 1, scaling); * *-----------------------------------------------------------------------* * int WebRtcIsac_Fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass, * size_t nSpan, int iSign, size_t max_factors, * size_t max_perm); * * RE, IM - see above documentation * * Although there is no limit on the number of dimensions, WebRtcIsac_Fftradix() must * be called once for each dimension, but the calls may be in any order. * * NTOTAL = the total number of complex data values * NPASS = the dimension of the current variable * NSPAN/NPASS = the spacing of consecutive data values while indexing the * current variable * ISIGN - see above documentation * * example: * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] * * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp); * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp); * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp); * * single-variate transform, * NTOTAL = N = NSPAN = (number of complex data values), * * WebRtcIsac_Fftradix (Re, Im, n, n, n, 1, maxf, maxp); * * The data can also be stored in a single array with alternating real and * imaginary parts, the magnitude of ISIGN is changed to 2 to give correct * indexing increment, and data [0] and data [1] used to pass the initial * addresses for the sequences of real and imaginary values, * * example: * REAL data [2*NTOTAL]; * WebRtcIsac_Fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp); * * for temporary allocation: * * MAX_FACTORS >= the maximum prime factor of NPASS * MAX_PERM >= the number of prime factors of NPASS. In addition, * if the square-free portion K of NPASS has two or more prime * factors, then MAX_PERM >= (K-1) * * storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS * has more than one square-free factor, the product of the square-free * factors must be <= 210 array storage for maximum prime factor of 23 the * following two constants should agree with the array dimensions. * *----------------------------------------------------------------------*/ #include #include #include "modules/third_party/fft/fft.h" /* double precision routine */ static int WebRtcIsac_Fftradix (double Re[], double Im[], size_t nTotal, size_t nPass, size_t nSpan, int isign, int max_factors, unsigned int max_perm, FFTstr *fftstate); #ifndef M_PI # define M_PI 3.14159265358979323846264338327950288 #endif #ifndef SIN60 # define SIN60 0.86602540378443865 /* sin(60 deg) */ # define COS72 0.30901699437494742 /* cos(72 deg) */ # define SIN72 0.95105651629515357 /* sin(72 deg) */ #endif # define REAL double # define FFTN WebRtcIsac_Fftn # define FFTNS "fftn" # define FFTRADIX WebRtcIsac_Fftradix # define FFTRADIXS "fftradix" int WebRtcIsac_Fftns(unsigned int ndim, const int dims[], double Re[], double Im[], int iSign, double scaling, FFTstr *fftstate) { size_t nSpan, nPass, nTotal; unsigned int i; int ret, max_factors, max_perm; /* * tally the number of elements in the data array * and determine the number of dimensions */ nTotal = 1; if (ndim && dims [0]) { for (i = 0; i < ndim; i++) { if (dims [i] <= 0) { return -1; } nTotal *= dims [i]; } } else { ndim = 0; for (i = 0; dims [i]; i++) { if (dims [i] <= 0) { return -1; } nTotal *= dims [i]; ndim++; } } /* determine maximum number of factors and permuations */ #if 1 /* * follow John Beale's example, just use the largest dimension and don't * worry about excess allocation. May be someone else will do it? */ max_factors = max_perm = 1; for (i = 0; i < ndim; i++) { nSpan = dims [i]; if ((int)nSpan > max_factors) { max_factors = (int)nSpan; } if ((int)nSpan > max_perm) { max_perm = (int)nSpan; } } #else /* use the constants used in the original Fortran code */ max_factors = 23; max_perm = 209; #endif /* loop over the dimensions: */ nPass = 1; for (i = 0; i < ndim; i++) { nSpan = dims [i]; nPass *= nSpan; ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign, max_factors, max_perm, fftstate); /* exit, clean-up already done */ if (ret) return ret; } /* Divide through by the normalizing constant: */ if (scaling && scaling != 1.0) { if (iSign < 0) iSign = -iSign; if (scaling < 0.0) { scaling = (double)nTotal; if (scaling < -1.0) scaling = sqrt (scaling); } scaling = 1.0 / scaling; /* multiply is often faster */ for (i = 0; i < nTotal; i += iSign) { Re [i] *= scaling; Im [i] *= scaling; } } return 0; } /* * singleton's mixed radix routine * * could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's * possible to make this a standalone function */ static int FFTRADIX (REAL Re[], REAL Im[], size_t nTotal, size_t nPass, size_t nSpan, int iSign, int max_factors, unsigned int max_perm, FFTstr *fftstate) { int ii, mfactor, kspan, ispan, inc; int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt; REAL radf; REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp; REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp; REAL *Rtmp = NULL; /* temp space for real part*/ REAL *Itmp = NULL; /* temp space for imaginary part */ REAL *Cos = NULL; /* Cosine values */ REAL *Sin = NULL; /* Sine values */ REAL s60 = SIN60; /* sin(60 deg) */ REAL c72 = COS72; /* cos(72 deg) */ REAL s72 = SIN72; /* sin(72 deg) */ REAL pi2 = M_PI; /* use PI first, 2 PI later */ fftstate->SpaceAlloced = 0; fftstate->MaxPermAlloced = 0; // initialize to avoid warnings k3 = c2 = c3 = s2 = s3 = 0.0; if (nPass < 2) return 0; /* allocate storage */ if (fftstate->SpaceAlloced < max_factors * sizeof (REAL)) { #ifdef SUN_BROKEN_REALLOC if (!fftstate->SpaceAlloced) /* first time */ { fftstate->SpaceAlloced = max_factors * sizeof (REAL); } else { #endif fftstate->SpaceAlloced = max_factors * sizeof (REAL); #ifdef SUN_BROKEN_REALLOC } #endif } else { /* allow full use of alloc'd space */ max_factors = fftstate->SpaceAlloced / sizeof (REAL); } if (fftstate->MaxPermAlloced < max_perm) { #ifdef SUN_BROKEN_REALLOC if (!fftstate->MaxPermAlloced) /* first time */ else #endif fftstate->MaxPermAlloced = max_perm; } else { /* allow full use of alloc'd space */ max_perm = fftstate->MaxPermAlloced; } /* assign pointers */ Rtmp = (REAL *) fftstate->Tmp0; Itmp = (REAL *) fftstate->Tmp1; Cos = (REAL *) fftstate->Tmp2; Sin = (REAL *) fftstate->Tmp3; /* * Function Body */ inc = iSign; if (iSign < 0) { s72 = -s72; s60 = -s60; pi2 = -pi2; inc = -inc; /* absolute value */ } /* adjust for strange increments */ nt = inc * (int)nTotal; ns = inc * (int)nSpan; kspan = ns; nn = nt - inc; jc = ns / (int)nPass; radf = pi2 * (double) jc; pi2 *= 2.0; /* use 2 PI from here on */ ii = 0; jf = 0; /* determine the factors of n */ mfactor = 0; k = (int)nPass; while (k % 16 == 0) { mfactor++; fftstate->factor [mfactor - 1] = 4; k /= 16; } j = 3; jj = 9; do { while (k % jj == 0) { mfactor++; fftstate->factor [mfactor - 1] = j; k /= jj; } j += 2; jj = j * j; } while (jj <= k); if (k <= 4) { kt = mfactor; fftstate->factor [mfactor] = k; if (k != 1) mfactor++; } else { if (k - (k / 4 << 2) == 0) { mfactor++; fftstate->factor [mfactor - 1] = 2; k /= 4; } kt = mfactor; j = 2; do { if (k % j == 0) { mfactor++; fftstate->factor [mfactor - 1] = j; k /= j; } j = ((j + 1) / 2 << 1) + 1; } while (j <= k); } if (kt) { j = kt; do { mfactor++; fftstate->factor [mfactor - 1] = fftstate->factor [j - 1]; j--; } while (j); } /* test that mfactors is in range */ if (mfactor > FFT_NFACTOR) { return -1; } /* compute fourier transform */ for (;;) { sd = radf / (double) kspan; cd = sin(sd); cd = 2.0 * cd * cd; sd = sin(sd + sd); kk = 0; ii++; switch (fftstate->factor [ii - 1]) { case 2: /* transform for factor of 2 (including rotation factor) */ kspan /= 2; k1 = kspan + 2; do { do { k2 = kk + kspan; ak = Re [k2]; bk = Im [k2]; Re [k2] = Re [kk] - ak; Im [k2] = Im [kk] - bk; Re [kk] += ak; Im [kk] += bk; kk = k2 + kspan; } while (kk < nn); kk -= nn; } while (kk < jc); if (kk >= kspan) goto Permute_Results_Label; /* exit infinite loop */ do { c1 = 1.0 - cd; s1 = sd; do { do { do { k2 = kk + kspan; ak = Re [kk] - Re [k2]; bk = Im [kk] - Im [k2]; Re [kk] += Re [k2]; Im [kk] += Im [k2]; Re [k2] = c1 * ak - s1 * bk; Im [k2] = s1 * ak + c1 * bk; kk = k2 + kspan; } while (kk < (nt-1)); k2 = kk - nt; c1 = -c1; kk = k1 - k2; } while (kk > k2); ak = c1 - (cd * c1 + sd * s1); s1 = sd * c1 - cd * s1 + s1; c1 = 2.0 - (ak * ak + s1 * s1); s1 *= c1; c1 *= ak; kk += jc; } while (kk < k2); k1 += inc + inc; kk = (k1 - kspan + 1) / 2 + jc - 1; } while (kk < (jc + jc)); break; case 4: /* transform for factor of 4 */ ispan = kspan; kspan /= 4; do { c1 = 1.0; s1 = 0.0; do { do { k1 = kk + kspan; k2 = k1 + kspan; k3 = k2 + kspan; akp = Re [kk] + Re [k2]; akm = Re [kk] - Re [k2]; ajp = Re [k1] + Re [k3]; ajm = Re [k1] - Re [k3]; bkp = Im [kk] + Im [k2]; bkm = Im [kk] - Im [k2]; bjp = Im [k1] + Im [k3]; bjm = Im [k1] - Im [k3]; Re [kk] = akp + ajp; Im [kk] = bkp + bjp; ajp = akp - ajp; bjp = bkp - bjp; if (iSign < 0) { akp = akm + bjm; bkp = bkm - ajm; akm -= bjm; bkm += ajm; } else { akp = akm - bjm; bkp = bkm + ajm; akm += bjm; bkm -= ajm; } /* avoid useless multiplies */ if (s1 == 0.0) { Re [k1] = akp; Re [k2] = ajp; Re [k3] = akm; Im [k1] = bkp; Im [k2] = bjp; Im [k3] = bkm; } else { Re [k1] = akp * c1 - bkp * s1; Re [k2] = ajp * c2 - bjp * s2; Re [k3] = akm * c3 - bkm * s3; Im [k1] = akp * s1 + bkp * c1; Im [k2] = ajp * s2 + bjp * c2; Im [k3] = akm * s3 + bkm * c3; } kk = k3 + kspan; } while (kk < nt); c2 = c1 - (cd * c1 + sd * s1); s1 = sd * c1 - cd * s1 + s1; c1 = 2.0 - (c2 * c2 + s1 * s1); s1 *= c1; c1 *= c2; /* values of c2, c3, s2, s3 that will get used next time */ c2 = c1 * c1 - s1 * s1; s2 = 2.0 * c1 * s1; c3 = c2 * c1 - s2 * s1; s3 = c2 * s1 + s2 * c1; kk = kk - nt + jc; } while (kk < kspan); kk = kk - kspan + inc; } while (kk < jc); if (kspan == jc) goto Permute_Results_Label; /* exit infinite loop */ break; default: /* transform for odd factors */ #ifdef FFT_RADIX4 return -1; break; #else /* FFT_RADIX4 */ k = fftstate->factor [ii - 1]; ispan = kspan; kspan /= k; switch (k) { case 3: /* transform for factor of 3 (optional code) */ do { do { k1 = kk + kspan; k2 = k1 + kspan; ak = Re [kk]; bk = Im [kk]; aj = Re [k1] + Re [k2]; bj = Im [k1] + Im [k2]; Re [kk] = ak + aj; Im [kk] = bk + bj; ak -= 0.5 * aj; bk -= 0.5 * bj; aj = (Re [k1] - Re [k2]) * s60; bj = (Im [k1] - Im [k2]) * s60; Re [k1] = ak - bj; Re [k2] = ak + bj; Im [k1] = bk + aj; Im [k2] = bk - aj; kk = k2 + kspan; } while (kk < (nn - 1)); kk -= nn; } while (kk < kspan); break; case 5: /* transform for factor of 5 (optional code) */ c2 = c72 * c72 - s72 * s72; s2 = 2.0 * c72 * s72; do { do { k1 = kk + kspan; k2 = k1 + kspan; k3 = k2 + kspan; k4 = k3 + kspan; akp = Re [k1] + Re [k4]; akm = Re [k1] - Re [k4]; bkp = Im [k1] + Im [k4]; bkm = Im [k1] - Im [k4]; ajp = Re [k2] + Re [k3]; ajm = Re [k2] - Re [k3]; bjp = Im [k2] + Im [k3]; bjm = Im [k2] - Im [k3]; aa = Re [kk]; bb = Im [kk]; Re [kk] = aa + akp + ajp; Im [kk] = bb + bkp + bjp; ak = akp * c72 + ajp * c2 + aa; bk = bkp * c72 + bjp * c2 + bb; aj = akm * s72 + ajm * s2; bj = bkm * s72 + bjm * s2; Re [k1] = ak - bj; Re [k4] = ak + bj; Im [k1] = bk + aj; Im [k4] = bk - aj; ak = akp * c2 + ajp * c72 + aa; bk = bkp * c2 + bjp * c72 + bb; aj = akm * s2 - ajm * s72; bj = bkm * s2 - bjm * s72; Re [k2] = ak - bj; Re [k3] = ak + bj; Im [k2] = bk + aj; Im [k3] = bk - aj; kk = k4 + kspan; } while (kk < (nn-1)); kk -= nn; } while (kk < kspan); break; default: if (k != jf) { jf = k; s1 = pi2 / (double) k; c1 = cos(s1); s1 = sin(s1); if (jf > max_factors){ return -1; } Cos [jf - 1] = 1.0; Sin [jf - 1] = 0.0; j = 1; do { Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1; Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1; k--; Cos [k - 1] = Cos [j - 1]; Sin [k - 1] = -Sin [j - 1]; j++; } while (j < k); } do { do { k1 = kk; k2 = kk + ispan; ak = aa = Re [kk]; bk = bb = Im [kk]; j = 1; k1 += kspan; do { k2 -= kspan; j++; Rtmp [j - 1] = Re [k1] + Re [k2]; ak += Rtmp [j - 1]; Itmp [j - 1] = Im [k1] + Im [k2]; bk += Itmp [j - 1]; j++; Rtmp [j - 1] = Re [k1] - Re [k2]; Itmp [j - 1] = Im [k1] - Im [k2]; k1 += kspan; } while (k1 < k2); Re [kk] = ak; Im [kk] = bk; k1 = kk; k2 = kk + ispan; j = 1; do { k1 += kspan; k2 -= kspan; jj = j; ak = aa; bk = bb; aj = 0.0; bj = 0.0; k = 1; do { k++; ak += Rtmp [k - 1] * Cos [jj - 1]; bk += Itmp [k - 1] * Cos [jj - 1]; k++; aj += Rtmp [k - 1] * Sin [jj - 1]; bj += Itmp [k - 1] * Sin [jj - 1]; jj += j; if (jj > jf) { jj -= jf; } } while (k < jf); k = jf - j; Re [k1] = ak - bj; Im [k1] = bk + aj; Re [k2] = ak + bj; Im [k2] = bk - aj; j++; } while (j < k); kk += ispan; } while (kk < nn); kk -= nn; } while (kk < kspan); break; } /* multiply by rotation factor (except for factors of 2 and 4) */ if (ii == mfactor) goto Permute_Results_Label; /* exit infinite loop */ kk = jc; do { c2 = 1.0 - cd; s1 = sd; do { c1 = c2; s2 = s1; kk += kspan; do { do { ak = Re [kk]; Re [kk] = c2 * ak - s2 * Im [kk]; Im [kk] = s2 * ak + c2 * Im [kk]; kk += ispan; } while (kk < nt); ak = s1 * s2; s2 = s1 * c2 + c1 * s2; c2 = c1 * c2 - ak; kk = kk - nt + kspan; } while (kk < ispan); c2 = c1 - (cd * c1 + sd * s1); s1 += sd * c1 - cd * s1; c1 = 2.0 - (c2 * c2 + s1 * s1); s1 *= c1; c2 *= c1; kk = kk - ispan + jc; } while (kk < kspan); kk = kk - kspan + jc + inc; } while (kk < (jc + jc)); break; #endif /* FFT_RADIX4 */ } } /* permute the results to normal order---done in two stages */ /* permutation for square factors of n */ Permute_Results_Label: fftstate->Perm [0] = ns; if (kt) { k = kt + kt + 1; if (mfactor < k) k--; j = 1; fftstate->Perm [k] = jc; do { fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1]; fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1]; j++; k--; } while (j < k); k3 = fftstate->Perm [k]; kspan = fftstate->Perm [1]; kk = jc; k2 = kspan; j = 1; if (nPass != nTotal) { /* permutation for multivariate transform */ Permute_Multi_Label: do { do { k = kk + jc; do { /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; kk += inc; k2 += inc; } while (kk < (k-1)); kk += ns - jc; k2 += ns - jc; } while (kk < (nt-1)); k2 = k2 - nt + kspan; kk = kk - nt + jc; } while (k2 < (ns-1)); do { do { k2 -= fftstate->Perm [j - 1]; j++; k2 = fftstate->Perm [j] + k2; } while (k2 > fftstate->Perm [j - 1]); j = 1; do { if (kk < (k2-1)) goto Permute_Multi_Label; kk += jc; k2 += kspan; } while (k2 < (ns-1)); } while (kk < (ns-1)); } else { /* permutation for single-variate transform (optional code) */ Permute_Single_Label: do { /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; kk += inc; k2 += kspan; } while (k2 < (ns-1)); do { do { k2 -= fftstate->Perm [j - 1]; j++; k2 = fftstate->Perm [j] + k2; } while (k2 >= fftstate->Perm [j - 1]); j = 1; do { if (kk < k2) goto Permute_Single_Label; kk += inc; k2 += kspan; } while (k2 < (ns-1)); } while (kk < (ns-1)); } jc = k3; } if ((kt << 1) + 1 >= mfactor) return 0; ispan = fftstate->Perm [kt]; /* permutation for square-free factors of n */ j = mfactor - kt; fftstate->factor [j] = 1; do { fftstate->factor [j - 1] *= fftstate->factor [j]; j--; } while (j != kt); kt++; nn = fftstate->factor [kt - 1] - 1; if (nn > (int) max_perm) { return -1; } j = jj = 0; for (;;) { k = kt + 1; k2 = fftstate->factor [kt - 1]; kk = fftstate->factor [k - 1]; j++; if (j > nn) break; /* exit infinite loop */ jj += kk; while (jj >= k2) { jj -= k2; k2 = kk; k++; kk = fftstate->factor [k - 1]; jj += kk; } fftstate->Perm [j - 1] = jj; } /* determine the permutation cycles of length greater than 1 */ j = 0; for (;;) { do { j++; kk = fftstate->Perm [j - 1]; } while (kk < 0); if (kk != j) { do { k = kk; kk = fftstate->Perm [k - 1]; fftstate->Perm [k - 1] = -kk; } while (kk != j); k3 = kk; } else { fftstate->Perm [j - 1] = -j; if (j == nn) break; /* exit infinite loop */ } } max_factors *= inc; /* reorder a and b, following the permutation cycles */ for (;;) { j = k3 + 1; nt -= ispan; ii = nt - inc + 1; if (nt < 0) break; /* exit infinite loop */ do { do { j--; } while (fftstate->Perm [j - 1] < 0); jj = jc; do { kspan = jj; if (jj > max_factors) { kspan = max_factors; } jj -= kspan; k = fftstate->Perm [j - 1]; kk = jc * k + ii + jj; k1 = kk + kspan - 1; k2 = 0; do { k2++; Rtmp [k2 - 1] = Re [k1]; Itmp [k2 - 1] = Im [k1]; k1 -= inc; } while (k1 != (kk-1)); do { k1 = kk + kspan - 1; k2 = k1 - jc * (k + fftstate->Perm [k - 1]); k = -fftstate->Perm [k - 1]; do { Re [k1] = Re [k2]; Im [k1] = Im [k2]; k1 -= inc; k2 -= inc; } while (k1 != (kk-1)); kk = k2 + 1; } while (k != j); k1 = kk + kspan - 1; k2 = 0; do { k2++; Re [k1] = Rtmp [k2 - 1]; Im [k1] = Itmp [k2 - 1]; k1 -= inc; } while (k1 != (kk-1)); } while (jj); } while (j != 1); } return 0; /* exit point here */ } /* ---------------------- end-of-file (c source) ---------------------- */