16f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin/* 26f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Copyright(c)1995,97 Mark Olesen <olesen@me.QueensU.CA> 36f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Queen's Univ at Kingston (Canada) 46f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 56f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Permission to use, copy, modify, and distribute this software for 66f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * any purpose without fee is hereby granted, provided that this 76f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * entire notice is included in all copies of any software which is 86f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * or includes a copy or modification of this software and in all 96f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * copies of the supporting documentation for such software. 106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR 126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S 136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY 146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS 156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * FITNESS FOR ANY PARTICULAR PURPOSE. 166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * All of which is to say that you can do what you like with this 186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * source code provided you don't try to sell it as your own and you 196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * include an unaltered copy of this message (including the 206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * copyright). 216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * It is also implicitly understood that bug fixes and improvements 236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * should make their way back to the general Internet community so 246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * that everyone benefits. 256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Changes: 276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Trivial type modifications by the WebRTC authors. 286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin */ 296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin/* 326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * File: 336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftn.c 346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Public: 366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftn / fftnf (); 376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Private: 396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftradix / fftradixf (); 406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Descript: 426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * multivariate complex Fourier transform, computed in place 436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * using mixed-radix Fast Fourier Transform algorithm. 446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Fortran code by: 466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * RC Singleton, Stanford Research Institute, Sept. 1968 476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * translated by f2c (version 19950721). 496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * int WebRtcIsac_Fftn (int ndim, const int dims[], REAL Re[], REAL Im[], 516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * int iSign, double scaling); 526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * NDIM = the total number dimensions 546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * DIMS = a vector of array sizes 556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * if NDIM is zero then DIMS must be zero-terminated 566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * RE and IM hold the real and imaginary components of the data, and return 586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * the resulting real and imaginary Fourier coefficients. Multidimensional 596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * data *must* be allocated contiguously. There is no limit on the number 606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * of dimensions. 616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * ISIGN = the sign of the complex exponential (ie, forward or inverse FFT) 636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * the magnitude of ISIGN (normally 1) is used to determine the 646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * correct indexing increment (see below). 656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * SCALING = normalizing constant by which the final result is *divided* 676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * if SCALING == -1, normalize by total dimension of the transform 686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * if SCALING < -1, normalize by the square-root of the total dimension 696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * example: 716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] 726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * int dims[3] = {n1,n2,n3} 746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftn (3, dims, Re, Im, 1, scaling); 756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin *-----------------------------------------------------------------------* 776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * int WebRtcIsac_Fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass, 786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * size_t nSpan, int iSign, size_t max_factors, 796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * size_t max_perm); 806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * RE, IM - see above documentation 826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Although there is no limit on the number of dimensions, WebRtcIsac_Fftradix() must 846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * be called once for each dimension, but the calls may be in any order. 856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * NTOTAL = the total number of complex data values 876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * NPASS = the dimension of the current variable 886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * NSPAN/NPASS = the spacing of consecutive data values while indexing the 896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * current variable 906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * ISIGN - see above documentation 916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * example: 936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] 946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp); 966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp); 976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp); 986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * single-variate transform, 1006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * NTOTAL = N = NSPAN = (number of complex data values), 1016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftradix (Re, Im, n, n, n, 1, maxf, maxp); 1036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * The data can also be stored in a single array with alternating real and 1056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * imaginary parts, the magnitude of ISIGN is changed to 2 to give correct 1066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * indexing increment, and data [0] and data [1] used to pass the initial 1076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * addresses for the sequences of real and imaginary values, 1086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * example: 1106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * REAL data [2*NTOTAL]; 1116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * WebRtcIsac_Fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp); 1126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * for temporary allocation: 1146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * MAX_FACTORS >= the maximum prime factor of NPASS 1166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * MAX_PERM >= the number of prime factors of NPASS. In addition, 1176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * if the square-free portion K of NPASS has two or more prime 1186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * factors, then MAX_PERM >= (K-1) 1196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS 1216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * has more than one square-free factor, the product of the square-free 1226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * factors must be <= 210 array storage for maximum prime factor of 23 the 1236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * following two constants should agree with the array dimensions. 1246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 1256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin *----------------------------------------------------------------------*/ 1266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#include "fft.h" 1276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#include <stdlib.h> 1296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#include <math.h> 1306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin/* double precision routine */ 1346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkinstatic int 1356f12fff925188ced26e518cd2252aff3e93bb04eAlexander GutkinWebRtcIsac_Fftradix (double Re[], double Im[], 1366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin size_t nTotal, size_t nPass, size_t nSpan, int isign, 1376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int max_factors, unsigned int max_perm, 1386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin FFTstr *fftstate); 1396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#ifndef M_PI 1436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define M_PI 3.14159265358979323846264338327950288 1446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif 1456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#ifndef SIN60 1476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define SIN60 0.86602540378443865 /* sin(60 deg) */ 1486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define COS72 0.30901699437494742 /* cos(72 deg) */ 1496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define SIN72 0.95105651629515357 /* sin(72 deg) */ 1506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif 1516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define REAL double 1536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define FFTN WebRtcIsac_Fftn 1546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define FFTNS "fftn" 1556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define FFTRADIX WebRtcIsac_Fftradix 1566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin# define FFTRADIXS "fftradix" 1576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkinint WebRtcIsac_Fftns(unsigned int ndim, const int dims[], 1606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin double Re[], 1616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin double Im[], 1626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int iSign, 1636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin double scaling, 1646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin FFTstr *fftstate) 1656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin{ 1666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin size_t nSpan, nPass, nTotal; 1686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin unsigned int i; 1696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int ret, max_factors, max_perm; 1706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 1716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* 1726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * tally the number of elements in the data array 1736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * and determine the number of dimensions 1746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin */ 1756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nTotal = 1; 1766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (ndim && dims [0]) 1776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 1786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (i = 0; i < ndim; i++) 1796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 1806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (dims [i] <= 0) 1816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 1826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 1836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 1846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nTotal *= dims [i]; 1856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 1866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 1876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin else 1886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 1896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ndim = 0; 1906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (i = 0; dims [i]; i++) 1916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 1926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (dims [i] <= 0) 1936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 1946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 1956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 1966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nTotal *= dims [i]; 1976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ndim++; 1986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 1996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* determine maximum number of factors and permuations */ 2026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#if 1 2036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* 2046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * follow John Beale's example, just use the largest dimension and don't 2056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * worry about excess allocation. May be someone else will do it? 2066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin */ 2076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_factors = max_perm = 1; 2086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (i = 0; i < ndim; i++) 2096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nSpan = dims [i]; 2116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if ((int)nSpan > max_factors) 2126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_factors = (int)nSpan; 2146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if ((int)nSpan > max_perm) 2166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_perm = (int)nSpan; 2186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#else 2216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* use the constants used in the original Fortran code */ 2226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_factors = 23; 2236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_perm = 209; 2246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif 2256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* loop over the dimensions: */ 2266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nPass = 1; 2276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (i = 0; i < ndim; i++) 2286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nSpan = dims [i]; 2306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nPass *= nSpan; 2316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign, 2326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_factors, max_perm, fftstate); 2336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* exit, clean-up already done */ 2346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (ret) 2356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return ret; 2366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* Divide through by the normalizing constant: */ 2396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (scaling && scaling != 1.0) 2406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (iSign < 0) iSign = -iSign; 2426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (scaling < 0.0) 2436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin scaling = (double)nTotal; 2456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (scaling < -1.0) 2466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin scaling = sqrt (scaling); 2476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin scaling = 1.0 / scaling; /* multiply is often faster */ 2496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (i = 0; i < nTotal; i += iSign) 2506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 2516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [i] *= scaling; 2526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [i] *= scaling; 2536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 2556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return 0; 2566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin} 2576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin/* 2596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * singleton's mixed radix routine 2606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * 2616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's 2626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * possible to make this a standalone function 2636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin */ 2646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkinstatic int FFTRADIX (REAL Re[], 2666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL Im[], 2676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin size_t nTotal, 2686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin size_t nPass, 2696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin size_t nSpan, 2706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int iSign, 2716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int max_factors, 2726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin unsigned int max_perm, 2736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin FFTstr *fftstate) 2746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin{ 2756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int ii, mfactor, kspan, ispan, inc; 2766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt; 2776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL radf; 2806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp; 2816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp; 2826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL *Rtmp = NULL; /* temp space for real part*/ 2846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL *Itmp = NULL; /* temp space for imaginary part */ 2856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL *Cos = NULL; /* Cosine values */ 2866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL *Sin = NULL; /* Sine values */ 2876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL s60 = SIN60; /* sin(60 deg) */ 2896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL c72 = COS72; /* cos(72 deg) */ 2906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL s72 = SIN72; /* sin(72 deg) */ 2916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin REAL pi2 = M_PI; /* use PI first, 2 PI later */ 2926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->SpaceAlloced = 0; 2956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->MaxPermAlloced = 0; 2966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 2986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin // initialize to avoid warnings 2996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k3 = c2 = c3 = s2 = s3 = 0.0; 3006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (nPass < 2) 3026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return 0; 3036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* allocate storage */ 3056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (fftstate->SpaceAlloced < max_factors * sizeof (REAL)) 3066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 3076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#ifdef SUN_BROKEN_REALLOC 3086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (!fftstate->SpaceAlloced) /* first time */ 3096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 3106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->SpaceAlloced = max_factors * sizeof (REAL); 3116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin else 3136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 3146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif 3156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->SpaceAlloced = max_factors * sizeof (REAL); 3166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#ifdef SUN_BROKEN_REALLOC 3176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif 3196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin else 3216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 3226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* allow full use of alloc'd space */ 3236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_factors = fftstate->SpaceAlloced / sizeof (REAL); 3246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (fftstate->MaxPermAlloced < max_perm) 3266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 3276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#ifdef SUN_BROKEN_REALLOC 3286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (!fftstate->MaxPermAlloced) /* first time */ 3296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin else 3306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif 3316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->MaxPermAlloced = max_perm; 3326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin else 3346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 3356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* allow full use of alloc'd space */ 3366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_perm = fftstate->MaxPermAlloced; 3376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (fftstate->Tmp0 == NULL || fftstate->Tmp1 == NULL || fftstate->Tmp2 == NULL || fftstate->Tmp3 == NULL 3396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin || fftstate->Perm == NULL) { 3406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 3416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* assign pointers */ 3446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Rtmp = (REAL *) fftstate->Tmp0; 3456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Itmp = (REAL *) fftstate->Tmp1; 3466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Cos = (REAL *) fftstate->Tmp2; 3476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Sin = (REAL *) fftstate->Tmp3; 3486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* 3506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin * Function Body 3516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin */ 3526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin inc = iSign; 3536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (iSign < 0) { 3546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s72 = -s72; 3556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s60 = -s60; 3566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin pi2 = -pi2; 3576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin inc = -inc; /* absolute value */ 3586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* adjust for strange increments */ 3616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nt = inc * (int)nTotal; 3626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ns = inc * (int)nSpan; 3636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan = ns; 3646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nn = nt - inc; 3666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jc = ns / (int)nPass; 3676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin radf = pi2 * (double) jc; 3686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin pi2 *= 2.0; /* use 2 PI from here on */ 3696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 3706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ii = 0; 3716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jf = 0; 3726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* determine the factors of n */ 3736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor = 0; 3746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = (int)nPass; 3756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin while (k % 16 == 0) { 3766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor++; 3776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [mfactor - 1] = 4; 3786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k /= 16; 3796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 3; 3816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj = 9; 3826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 3836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin while (k % jj == 0) { 3846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor++; 3856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [mfactor - 1] = j; 3866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k /= jj; 3876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 3886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j += 2; 3896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj = j * j; 3906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (jj <= k); 3916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (k <= 4) { 3926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kt = mfactor; 3936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [mfactor] = k; 3946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (k != 1) 3956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor++; 3966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } else { 3976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (k - (k / 4 << 2) == 0) { 3986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor++; 3996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [mfactor - 1] = 2; 4006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k /= 4; 4016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 4026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kt = mfactor; 4036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 2; 4046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (k % j == 0) { 4066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor++; 4076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [mfactor - 1] = j; 4086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k /= j; 4096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 4106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = ((j + 1) / 2 << 1) + 1; 4116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j <= k); 4126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 4136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kt) { 4146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = kt; 4156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin mfactor++; 4176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [mfactor - 1] = fftstate->factor [j - 1]; 4186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j--; 4196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j); 4206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 4216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 4226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* test that mfactors is in range */ 4236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (mfactor > NFACTOR) 4246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin { 4256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 4266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 4276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 4286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* compute fourier transform */ 4296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (;;) { 4306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin sd = radf / (double) kspan; 4316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin cd = sin(sd); 4326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin cd = 2.0 * cd * cd; 4336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin sd = sin(sd + sd); 4346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = 0; 4356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ii++; 4366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 4376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin switch (fftstate->factor [ii - 1]) { 4386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin case 2: 4396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* transform for factor of 2 (including rotation factor) */ 4406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan /= 2; 4416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kspan + 2; 4426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kk + kspan; 4456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = Re [k2]; 4466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = Im [k2]; 4476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = Re [kk] - ak; 4486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = Im [kk] - bk; 4496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] += ak; 4506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] += bk; 4516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k2 + kspan; 4526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < nn); 4536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk -= nn; 4546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < jc); 4556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kk >= kspan) 4566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin goto Permute_Results_Label; /* exit infinite loop */ 4576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = 1.0 - cd; 4596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = sd; 4606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kk + kspan; 4646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = Re [kk] - Re [k2]; 4656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = Im [kk] - Im [k2]; 4666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] += Re [k2]; 4676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] += Im [k2]; 4686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = c1 * ak - s1 * bk; 4696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = s1 * ak + c1 * bk; 4706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k2 + kspan; 4716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (nt-1)); 4726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kk - nt; 4736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = -c1; 4746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k1 - k2; 4756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk > k2); 4766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = c1 - (cd * c1 + sd * s1); 4776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = sd * c1 - cd * s1 + s1; 4786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = 2.0 - (ak * ak + s1 * s1); 4796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 *= c1; 4806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 *= ak; 4816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += jc; 4826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < k2); 4836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 += inc + inc; 4846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = (k1 - kspan + 1) / 2 + jc - 1; 4856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (jc + jc)); 4866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 4876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 4886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin case 4: /* transform for factor of 4 */ 4896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ispan = kspan; 4906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan /= 4; 4916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 4926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = 1.0; 4946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = 0.0; 4956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 4976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk + kspan; 4986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = k1 + kspan; 4996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k3 = k2 + kspan; 5006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akp = Re [kk] + Re [k2]; 5016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akm = Re [kk] - Re [k2]; 5026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ajp = Re [k1] + Re [k3]; 5036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ajm = Re [k1] - Re [k3]; 5046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkp = Im [kk] + Im [k2]; 5056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkm = Im [kk] - Im [k2]; 5066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bjp = Im [k1] + Im [k3]; 5076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bjm = Im [k1] - Im [k3]; 5086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] = akp + ajp; 5096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] = bkp + bjp; 5106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ajp = akp - ajp; 5116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bjp = bkp - bjp; 5126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (iSign < 0) { 5136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akp = akm + bjm; 5146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkp = bkm - ajm; 5156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akm -= bjm; 5166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkm += ajm; 5176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } else { 5186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akp = akm - bjm; 5196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkp = bkm + ajm; 5206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akm += bjm; 5216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkm -= ajm; 5226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 5236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* avoid useless multiplies */ 5246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (s1 == 0.0) { 5256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = akp; 5266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = ajp; 5276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k3] = akm; 5286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = bkp; 5296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = bjp; 5306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k3] = bkm; 5316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } else { 5326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = akp * c1 - bkp * s1; 5336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = ajp * c2 - bjp * s2; 5346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k3] = akm * c3 - bkm * s3; 5356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = akp * s1 + bkp * c1; 5366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = ajp * s2 + bjp * c2; 5376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k3] = akm * s3 + bkm * c3; 5386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 5396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k3 + kspan; 5406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < nt); 5416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 5426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 = c1 - (cd * c1 + sd * s1); 5436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = sd * c1 - cd * s1 + s1; 5446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = 2.0 - (c2 * c2 + s1 * s1); 5456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 *= c1; 5466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 *= c2; 5476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* values of c2, c3, s2, s3 that will get used next time */ 5486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 = c1 * c1 - s1 * s1; 5496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s2 = 2.0 * c1 * s1; 5506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c3 = c2 * c1 - s2 * s1; 5516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s3 = c2 * s1 + s2 * c1; 5526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = kk - nt + jc; 5536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < kspan); 5546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = kk - kspan + inc; 5556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < jc); 5566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kspan == jc) 5576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin goto Permute_Results_Label; /* exit infinite loop */ 5586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 5596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 5606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin default: 5616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* transform for odd factors */ 5626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#ifdef FFT_RADIX4 5636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 5646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 5656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#else /* FFT_RADIX4 */ 5666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = fftstate->factor [ii - 1]; 5676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ispan = kspan; 5686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan /= k; 5696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 5706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin switch (k) { 5716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin case 3: /* transform for factor of 3 (optional code) */ 5726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 5736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 5746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk + kspan; 5756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = k1 + kspan; 5766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = Re [kk]; 5776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = Im [kk]; 5786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aj = Re [k1] + Re [k2]; 5796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bj = Im [k1] + Im [k2]; 5806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] = ak + aj; 5816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] = bk + bj; 5826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak -= 0.5 * aj; 5836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk -= 0.5 * bj; 5846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aj = (Re [k1] - Re [k2]) * s60; 5856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bj = (Im [k1] - Im [k2]) * s60; 5866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = ak - bj; 5876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = ak + bj; 5886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = bk + aj; 5896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = bk - aj; 5906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k2 + kspan; 5916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (nn - 1)); 5926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk -= nn; 5936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < kspan); 5946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 5956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 5966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin case 5: /* transform for factor of 5 (optional code) */ 5976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 = c72 * c72 - s72 * s72; 5986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s2 = 2.0 * c72 * s72; 5996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk + kspan; 6026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = k1 + kspan; 6036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k3 = k2 + kspan; 6046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k4 = k3 + kspan; 6056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akp = Re [k1] + Re [k4]; 6066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin akm = Re [k1] - Re [k4]; 6076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkp = Im [k1] + Im [k4]; 6086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bkm = Im [k1] - Im [k4]; 6096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ajp = Re [k2] + Re [k3]; 6106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ajm = Re [k2] - Re [k3]; 6116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bjp = Im [k2] + Im [k3]; 6126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bjm = Im [k2] - Im [k3]; 6136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aa = Re [kk]; 6146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bb = Im [kk]; 6156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] = aa + akp + ajp; 6166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] = bb + bkp + bjp; 6176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = akp * c72 + ajp * c2 + aa; 6186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = bkp * c72 + bjp * c2 + bb; 6196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aj = akm * s72 + ajm * s2; 6206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bj = bkm * s72 + bjm * s2; 6216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = ak - bj; 6226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k4] = ak + bj; 6236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = bk + aj; 6246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k4] = bk - aj; 6256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = akp * c2 + ajp * c72 + aa; 6266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = bkp * c2 + bjp * c72 + bb; 6276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aj = akm * s2 - ajm * s72; 6286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bj = bkm * s2 - bjm * s72; 6296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = ak - bj; 6306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k3] = ak + bj; 6316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = bk + aj; 6326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k3] = bk - aj; 6336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k4 + kspan; 6346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (nn-1)); 6356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk -= nn; 6366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < kspan); 6376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 6386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 6396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin default: 6406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (k != jf) { 6416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jf = k; 6426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = pi2 / (double) k; 6436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = cos(s1); 6446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = sin(s1); 6456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (jf > max_factors){ 6466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 6476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 6486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Cos [jf - 1] = 1.0; 6496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Sin [jf - 1] = 0.0; 6506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 6516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1; 6536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1; 6546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k--; 6556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Cos [k - 1] = Cos [j - 1]; 6566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Sin [k - 1] = -Sin [j - 1]; 6576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 6586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j < k); 6596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 6606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk; 6636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kk + ispan; 6646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = aa = Re [kk]; 6656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = bb = Im [kk]; 6666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 6676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 += kspan; 6686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 -= kspan; 6706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 6716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Rtmp [j - 1] = Re [k1] + Re [k2]; 6726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak += Rtmp [j - 1]; 6736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Itmp [j - 1] = Im [k1] + Im [k2]; 6746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk += Itmp [j - 1]; 6756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 6766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Rtmp [j - 1] = Re [k1] - Re [k2]; 6776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Itmp [j - 1] = Im [k1] - Im [k2]; 6786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 += kspan; 6796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k1 < k2); 6806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] = ak; 6816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] = bk; 6826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk; 6836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kk + ispan; 6846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 6856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 += kspan; 6876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 -= kspan; 6886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj = j; 6896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = aa; 6906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = bb; 6916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aj = 0.0; 6926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bj = 0.0; 6936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = 1; 6946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 6956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k++; 6966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak += Rtmp [k - 1] * Cos [jj - 1]; 6976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk += Itmp [k - 1] * Cos [jj - 1]; 6986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k++; 6996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin aj += Rtmp [k - 1] * Sin [jj - 1]; 7006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bj += Itmp [k - 1] * Sin [jj - 1]; 7016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj += j; 7026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (jj > jf) { 7036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj -= jf; 7046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 7056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k < jf); 7066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = jf - j; 7076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = ak - bj; 7086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = bk + aj; 7096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k2] = ak + bj; 7106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k2] = bk - aj; 7116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 7126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j < k); 7136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += ispan; 7146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < nn); 7156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk -= nn; 7166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < kspan); 7176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 7186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 7196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 7206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* multiply by rotation factor (except for factors of 2 and 4) */ 7216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (ii == mfactor) 7226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin goto Permute_Results_Label; /* exit infinite loop */ 7236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = jc; 7246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 = 1.0 - cd; 7266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 = sd; 7276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = c2; 7296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s2 = s1; 7306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += kspan; 7316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = Re [kk]; 7346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [kk] = c2 * ak - s2 * Im [kk]; 7356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [kk] = s2 * ak + c2 * Im [kk]; 7366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += ispan; 7376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < nt); 7386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = s1 * s2; 7396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s2 = s1 * c2 + c1 * s2; 7406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 = c1 * c2 - ak; 7416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = kk - nt + kspan; 7426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < ispan); 7436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 = c1 - (cd * c1 + sd * s1); 7446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 += sd * c1 - cd * s1; 7456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c1 = 2.0 - (c2 * c2 + s1 * s1); 7466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin s1 *= c1; 7476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin c2 *= c1; 7486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = kk - ispan + jc; 7496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < kspan); 7506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = kk - kspan + jc + inc; 7516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (jc + jc)); 7526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; 7536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin#endif /* FFT_RADIX4 */ 7546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 7556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 7566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 7576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* permute the results to normal order---done in two stages */ 7586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* permutation for square factors of n */ 7596f12fff925188ced26e518cd2252aff3e93bb04eAlexander GutkinPermute_Results_Label: 7606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [0] = ns; 7616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kt) { 7626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = kt + kt + 1; 7636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (mfactor < k) 7646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k--; 7656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 7666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [k] = jc; 7676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1]; 7696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1]; 7706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 7716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k--; 7726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j < k); 7736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k3 = fftstate->Perm [k]; 7746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan = fftstate->Perm [1]; 7756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = jc; 7766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kspan; 7776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 7786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (nPass != nTotal) { 7796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* permutation for multivariate transform */ 7806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Permute_Multi_Label: 7816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = kk + jc; 7846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ 7866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; 7876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; 7886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += inc; 7896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 += inc; 7906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (k-1)); 7916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += ns - jc; 7926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 += ns - jc; 7936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (nt-1)); 7946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = k2 - nt + kspan; 7956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = kk - nt + jc; 7966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k2 < (ns-1)); 7976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 7996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 -= fftstate->Perm [j - 1]; 8006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 8016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = fftstate->Perm [j] + k2; 8026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k2 > fftstate->Perm [j - 1]); 8036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 8046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kk < (k2-1)) 8066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin goto Permute_Multi_Label; 8076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += jc; 8086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 += kspan; 8096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k2 < (ns-1)); 8106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (ns-1)); 8116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } else { 8126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* permutation for single-variate transform (optional code) */ 8136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Permute_Single_Label: 8146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ 8166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; 8176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; 8186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += inc; 8196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 += kspan; 8206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k2 < (ns-1)); 8216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 -= fftstate->Perm [j - 1]; 8246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 8256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = fftstate->Perm [j] + k2; 8266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k2 >= fftstate->Perm [j - 1]); 8276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 1; 8286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kk < k2) 8306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin goto Permute_Single_Label; 8316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk += inc; 8326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 += kspan; 8336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k2 < (ns-1)); 8346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < (ns-1)); 8356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jc = k3; 8376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 8396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if ((kt << 1) + 1 >= mfactor) 8406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return 0; 8416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ispan = fftstate->Perm [kt]; 8426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* permutation for square-free factors of n */ 8436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = mfactor - kt; 8446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [j] = 1; 8456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->factor [j - 1] *= fftstate->factor [j]; 8476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j--; 8486f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j != kt); 8496f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kt++; 8506f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nn = fftstate->factor [kt - 1] - 1; 8516f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (nn > (int) max_perm) { 8526f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return -1; 8536f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8546f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = jj = 0; 8556f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (;;) { 8566f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = kt + 1; 8576f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = fftstate->factor [kt - 1]; 8586f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = fftstate->factor [k - 1]; 8596f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 8606f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (j > nn) 8616f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; /* exit infinite loop */ 8626f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj += kk; 8636f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin while (jj >= k2) { 8646f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj -= k2; 8656f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = kk; 8666f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k++; 8676f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = fftstate->factor [k - 1]; 8686f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj += kk; 8696f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8706f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [j - 1] = jj; 8716f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8726f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* determine the permutation cycles of length greater than 1 */ 8736f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = 0; 8746f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (;;) { 8756f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8766f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j++; 8776f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = fftstate->Perm [j - 1]; 8786f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk < 0); 8796f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (kk != j) { 8806f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 8816f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = kk; 8826f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = fftstate->Perm [k - 1]; 8836f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [k - 1] = -kk; 8846f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (kk != j); 8856f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k3 = kk; 8866f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } else { 8876f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin fftstate->Perm [j - 1] = -j; 8886f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (j == nn) 8896f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; /* exit infinite loop */ 8906f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8916f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 8926f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin max_factors *= inc; 8936f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin /* reorder a and b, following the permutation cycles */ 8946f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin for (;;) { 8956f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j = k3 + 1; 8966f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin nt -= ispan; 8976f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin ii = nt - inc + 1; 8986f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (nt < 0) 8996f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin break; /* exit infinite loop */ 9006f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9016f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9026f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin j--; 9036f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (fftstate->Perm [j - 1] < 0); 9046f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj = jc; 9056f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9066f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan = jj; 9076f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin if (jj > max_factors) { 9086f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kspan = max_factors; 9096f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 9106f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin jj -= kspan; 9116f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = fftstate->Perm [j - 1]; 9126f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = jc * k + ii + jj; 9136f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk + kspan - 1; 9146f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = 0; 9156f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9166f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2++; 9176f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Rtmp [k2 - 1] = Re [k1]; 9186f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Itmp [k2 - 1] = Im [k1]; 9196f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 -= inc; 9206f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k1 != (kk-1)); 9216f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9226f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk + kspan - 1; 9236f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = k1 - jc * (k + fftstate->Perm [k - 1]); 9246f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k = -fftstate->Perm [k - 1]; 9256f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9266f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = Re [k2]; 9276f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = Im [k2]; 9286f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 -= inc; 9296f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 -= inc; 9306f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k1 != (kk-1)); 9316f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin kk = k2 + 1; 9326f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k != j); 9336f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 = kk + kspan - 1; 9346f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2 = 0; 9356f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin do { 9366f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k2++; 9376f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Re [k1] = Rtmp [k2 - 1]; 9386f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin Im [k1] = Itmp [k2 - 1]; 9396f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin k1 -= inc; 9406f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (k1 != (kk-1)); 9416f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (jj); 9426f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } while (j != 1); 9436f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin } 9446f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin return 0; /* exit point here */ 9456f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin} 9466f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin/* ---------------------- end-of-file (c source) ---------------------- */ 9476f12fff925188ced26e518cd2252aff3e93bb04eAlexander Gutkin 948