1b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org/* 2b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Copyright(c)1995,97 Mark Olesen <olesen@me.QueensU.CA> 3b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Queen's Univ at Kingston (Canada) 4b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 5b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Permission to use, copy, modify, and distribute this software for 6b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * any purpose without fee is hereby granted, provided that this 7b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * entire notice is included in all copies of any software which is 8b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * or includes a copy or modification of this software and in all 9b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * copies of the supporting documentation for such software. 10b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 11b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR 12b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S 13b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY 14b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS 15b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * FITNESS FOR ANY PARTICULAR PURPOSE. 16b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 17b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * All of which is to say that you can do what you like with this 18b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * source code provided you don't try to sell it as your own and you 19b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * include an unaltered copy of this message (including the 20b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * copyright). 21b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 22b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * It is also implicitly understood that bug fixes and improvements 23b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * should make their way back to the general Internet community so 24b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * that everyone benefits. 25b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 26b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Changes: 27b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Trivial type modifications by the WebRTC authors. 28b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org */ 29b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 30b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 31b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org/* 32b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * File: 33b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftn.c 34b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 35b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Public: 36b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftn / fftnf (); 37b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 38b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Private: 39b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftradix / fftradixf (); 40b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 41b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Descript: 42b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * multivariate complex Fourier transform, computed in place 43b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * using mixed-radix Fast Fourier Transform algorithm. 44b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 45b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Fortran code by: 46b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * RC Singleton, Stanford Research Institute, Sept. 1968 47b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 48b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * translated by f2c (version 19950721). 49b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 50b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * int WebRtcIsac_Fftn (int ndim, const int dims[], REAL Re[], REAL Im[], 51b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * int iSign, double scaling); 52b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 53b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * NDIM = the total number dimensions 54b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * DIMS = a vector of array sizes 55b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * if NDIM is zero then DIMS must be zero-terminated 56b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 57b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * RE and IM hold the real and imaginary components of the data, and return 58b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * the resulting real and imaginary Fourier coefficients. Multidimensional 59b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * data *must* be allocated contiguously. There is no limit on the number 60b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * of dimensions. 61b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 62b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * ISIGN = the sign of the complex exponential (ie, forward or inverse FFT) 63b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * the magnitude of ISIGN (normally 1) is used to determine the 64b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * correct indexing increment (see below). 65b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 66b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * SCALING = normalizing constant by which the final result is *divided* 67b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * if SCALING == -1, normalize by total dimension of the transform 68b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * if SCALING < -1, normalize by the square-root of the total dimension 69b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 70b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * example: 71b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] 72b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 73b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * int dims[3] = {n1,n2,n3} 74b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftn (3, dims, Re, Im, 1, scaling); 75b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 76b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org *-----------------------------------------------------------------------* 77b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * int WebRtcIsac_Fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass, 78b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * size_t nSpan, int iSign, size_t max_factors, 79b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * size_t max_perm); 80b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 81b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * RE, IM - see above documentation 82b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 83b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Although there is no limit on the number of dimensions, WebRtcIsac_Fftradix() must 84b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * be called once for each dimension, but the calls may be in any order. 85b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 86b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * NTOTAL = the total number of complex data values 87b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * NPASS = the dimension of the current variable 88b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * NSPAN/NPASS = the spacing of consecutive data values while indexing the 89b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * current variable 90b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * ISIGN - see above documentation 91b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 92b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * example: 93b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] 94b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 95b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp); 96b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp); 97b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp); 98b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 99b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * single-variate transform, 100b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * NTOTAL = N = NSPAN = (number of complex data values), 101b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 102b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftradix (Re, Im, n, n, n, 1, maxf, maxp); 103b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 104b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * The data can also be stored in a single array with alternating real and 105b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * imaginary parts, the magnitude of ISIGN is changed to 2 to give correct 106b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * indexing increment, and data [0] and data [1] used to pass the initial 107b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * addresses for the sequences of real and imaginary values, 108b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 109b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * example: 110b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * REAL data [2*NTOTAL]; 111b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * WebRtcIsac_Fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp); 112b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 113b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * for temporary allocation: 114b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 115b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * MAX_FACTORS >= the maximum prime factor of NPASS 116b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * MAX_PERM >= the number of prime factors of NPASS. In addition, 117b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * if the square-free portion K of NPASS has two or more prime 118b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * factors, then MAX_PERM >= (K-1) 119b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 120b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS 121b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * has more than one square-free factor, the product of the square-free 122b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * factors must be <= 210 array storage for maximum prime factor of 23 the 123b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * following two constants should agree with the array dimensions. 124b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 125b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org *----------------------------------------------------------------------*/ 126b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#include "fft.h" 127b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 128b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#include <stdlib.h> 129b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#include <math.h> 130b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 131b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 132b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 133b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org/* double precision routine */ 134b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.orgstatic int 135b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.orgWebRtcIsac_Fftradix (double Re[], double Im[], 136b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org size_t nTotal, size_t nPass, size_t nSpan, int isign, 137b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int max_factors, unsigned int max_perm, 138b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org FFTstr *fftstate); 139b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 140b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 141b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 142b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#ifndef M_PI 143b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define M_PI 3.14159265358979323846264338327950288 144b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif 145b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 146b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#ifndef SIN60 147b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define SIN60 0.86602540378443865 /* sin(60 deg) */ 148b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define COS72 0.30901699437494742 /* cos(72 deg) */ 149b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define SIN72 0.95105651629515357 /* sin(72 deg) */ 150b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif 151b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 152b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define REAL double 153b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define FFTN WebRtcIsac_Fftn 154b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define FFTNS "fftn" 155b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define FFTRADIX WebRtcIsac_Fftradix 156b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org# define FFTRADIXS "fftradix" 157b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 158b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 159b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.orgint WebRtcIsac_Fftns(unsigned int ndim, const int dims[], 160b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org double Re[], 161b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org double Im[], 162b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int iSign, 163b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org double scaling, 164b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org FFTstr *fftstate) 165b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org{ 166b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 167b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org size_t nSpan, nPass, nTotal; 168b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org unsigned int i; 169b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int ret, max_factors, max_perm; 170b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 171b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* 172b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * tally the number of elements in the data array 173b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * and determine the number of dimensions 174b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org */ 175b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nTotal = 1; 176b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (ndim && dims [0]) 177b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 178b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (i = 0; i < ndim; i++) 179b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 180b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (dims [i] <= 0) 181b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 182b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return -1; 183b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 184b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nTotal *= dims [i]; 185b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 186b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 187b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org else 188b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 189b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ndim = 0; 190b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (i = 0; dims [i]; i++) 191b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 192b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (dims [i] <= 0) 193b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 194b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return -1; 195b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 196b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nTotal *= dims [i]; 197b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ndim++; 198b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 199b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 200b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 201b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* determine maximum number of factors and permuations */ 202b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#if 1 203b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* 204b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * follow John Beale's example, just use the largest dimension and don't 205b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * worry about excess allocation. May be someone else will do it? 206b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org */ 207b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_factors = max_perm = 1; 208b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (i = 0; i < ndim; i++) 209b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 210b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nSpan = dims [i]; 211b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if ((int)nSpan > max_factors) 212b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 213b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_factors = (int)nSpan; 214b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 215b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if ((int)nSpan > max_perm) 216b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 217b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_perm = (int)nSpan; 218b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 219b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 220b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#else 221b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* use the constants used in the original Fortran code */ 222b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_factors = 23; 223b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_perm = 209; 224b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif 225b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* loop over the dimensions: */ 226b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nPass = 1; 227b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (i = 0; i < ndim; i++) 228b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 229b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nSpan = dims [i]; 230b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nPass *= nSpan; 231b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign, 232b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_factors, max_perm, fftstate); 233b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* exit, clean-up already done */ 234b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (ret) 235b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return ret; 236b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 237b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 238b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* Divide through by the normalizing constant: */ 239b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (scaling && scaling != 1.0) 240b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 241b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (iSign < 0) iSign = -iSign; 242b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (scaling < 0.0) 243b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 244b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org scaling = (double)nTotal; 245b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (scaling < -1.0) 246b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org scaling = sqrt (scaling); 247b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 248b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org scaling = 1.0 / scaling; /* multiply is often faster */ 249b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (i = 0; i < nTotal; i += iSign) 250b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 251b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [i] *= scaling; 252b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [i] *= scaling; 253b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 254b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 255b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return 0; 256b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org} 257b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 258b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org/* 259b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * singleton's mixed radix routine 260b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * 261b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's 262b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * possible to make this a standalone function 263b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org */ 264b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 265b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.orgstatic int FFTRADIX (REAL Re[], 266b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL Im[], 267b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org size_t nTotal, 268b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org size_t nPass, 269b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org size_t nSpan, 270b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int iSign, 271b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int max_factors, 272b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org unsigned int max_perm, 273b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org FFTstr *fftstate) 274b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org{ 275b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int ii, mfactor, kspan, ispan, inc; 276b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt; 277b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 278b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 279b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL radf; 280b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp; 281b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp; 282b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 283b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL *Rtmp = NULL; /* temp space for real part*/ 284b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL *Itmp = NULL; /* temp space for imaginary part */ 285b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL *Cos = NULL; /* Cosine values */ 286b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL *Sin = NULL; /* Sine values */ 287b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 288b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL s60 = SIN60; /* sin(60 deg) */ 289b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL c72 = COS72; /* cos(72 deg) */ 290b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL s72 = SIN72; /* sin(72 deg) */ 291b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org REAL pi2 = M_PI; /* use PI first, 2 PI later */ 292b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 293b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 294b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->SpaceAlloced = 0; 295b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->MaxPermAlloced = 0; 296b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 297b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 298b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org // initialize to avoid warnings 299b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k3 = c2 = c3 = s2 = s3 = 0.0; 300b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 301b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (nPass < 2) 302b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return 0; 303b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 304b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* allocate storage */ 305b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (fftstate->SpaceAlloced < max_factors * sizeof (REAL)) 306b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 307b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#ifdef SUN_BROKEN_REALLOC 308b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (!fftstate->SpaceAlloced) /* first time */ 309b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 310b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->SpaceAlloced = max_factors * sizeof (REAL); 311b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 312b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org else 313b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 314b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif 315b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->SpaceAlloced = max_factors * sizeof (REAL); 316b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#ifdef SUN_BROKEN_REALLOC 317b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 318b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif 319b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 320b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org else 321b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 322b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* allow full use of alloc'd space */ 323b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_factors = fftstate->SpaceAlloced / sizeof (REAL); 324b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 325b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (fftstate->MaxPermAlloced < max_perm) 326b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 327b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#ifdef SUN_BROKEN_REALLOC 328b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (!fftstate->MaxPermAlloced) /* first time */ 329b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org else 330b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif 331b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->MaxPermAlloced = max_perm; 332b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 333b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org else 334b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 335b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* allow full use of alloc'd space */ 336b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_perm = fftstate->MaxPermAlloced; 337b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 338b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 339b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* assign pointers */ 340b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Rtmp = (REAL *) fftstate->Tmp0; 341b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Itmp = (REAL *) fftstate->Tmp1; 342b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Cos = (REAL *) fftstate->Tmp2; 343b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Sin = (REAL *) fftstate->Tmp3; 344b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 345b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* 346b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org * Function Body 347b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org */ 348b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org inc = iSign; 349b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (iSign < 0) { 350b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s72 = -s72; 351b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s60 = -s60; 352b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org pi2 = -pi2; 353b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org inc = -inc; /* absolute value */ 354b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 355b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 356b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* adjust for strange increments */ 357b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nt = inc * (int)nTotal; 358b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ns = inc * (int)nSpan; 359b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan = ns; 360b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 361b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nn = nt - inc; 362b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jc = ns / (int)nPass; 363b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org radf = pi2 * (double) jc; 364b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org pi2 *= 2.0; /* use 2 PI from here on */ 365b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 366b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ii = 0; 367b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jf = 0; 368b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* determine the factors of n */ 369b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor = 0; 370b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = (int)nPass; 371b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org while (k % 16 == 0) { 372b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor++; 373b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [mfactor - 1] = 4; 374b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k /= 16; 375b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 376b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 3; 377b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj = 9; 378b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 379b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org while (k % jj == 0) { 380b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor++; 381b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [mfactor - 1] = j; 382b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k /= jj; 383b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 384b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j += 2; 385b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj = j * j; 386b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (jj <= k); 387b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (k <= 4) { 388b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kt = mfactor; 389b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [mfactor] = k; 390b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (k != 1) 391b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor++; 392b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } else { 393b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (k - (k / 4 << 2) == 0) { 394b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor++; 395b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [mfactor - 1] = 2; 396b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k /= 4; 397b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 398b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kt = mfactor; 399b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 2; 400b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 401b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (k % j == 0) { 402b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor++; 403b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [mfactor - 1] = j; 404b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k /= j; 405b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 406b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = ((j + 1) / 2 << 1) + 1; 407b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j <= k); 408b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 409b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kt) { 410b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = kt; 411b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 412b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org mfactor++; 413b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [mfactor - 1] = fftstate->factor [j - 1]; 414b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j--; 415b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j); 416b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 417b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 418b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* test that mfactors is in range */ 419b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (mfactor > NFACTOR) 420b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org { 421b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return -1; 422b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 423b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 424b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* compute fourier transform */ 425b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (;;) { 426b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org sd = radf / (double) kspan; 427b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org cd = sin(sd); 428b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org cd = 2.0 * cd * cd; 429b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org sd = sin(sd + sd); 430b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = 0; 431b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ii++; 432b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 433b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org switch (fftstate->factor [ii - 1]) { 434b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org case 2: 435b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* transform for factor of 2 (including rotation factor) */ 436b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan /= 2; 437b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kspan + 2; 438b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 439b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 440b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kk + kspan; 441b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = Re [k2]; 442b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = Im [k2]; 443b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = Re [kk] - ak; 444b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = Im [kk] - bk; 445b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] += ak; 446b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] += bk; 447b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k2 + kspan; 448b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < nn); 449b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk -= nn; 450b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < jc); 451b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kk >= kspan) 452b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org goto Permute_Results_Label; /* exit infinite loop */ 453b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 454b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = 1.0 - cd; 455b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = sd; 456b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 457b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 458b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 459b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kk + kspan; 460b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = Re [kk] - Re [k2]; 461b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = Im [kk] - Im [k2]; 462b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] += Re [k2]; 463b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] += Im [k2]; 464b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = c1 * ak - s1 * bk; 465b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = s1 * ak + c1 * bk; 466b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k2 + kspan; 467b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (nt-1)); 468b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kk - nt; 469b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = -c1; 470b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k1 - k2; 471b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk > k2); 472b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = c1 - (cd * c1 + sd * s1); 473b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = sd * c1 - cd * s1 + s1; 474b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = 2.0 - (ak * ak + s1 * s1); 475b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 *= c1; 476b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 *= ak; 477b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += jc; 478b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < k2); 479b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 += inc + inc; 480b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = (k1 - kspan + 1) / 2 + jc - 1; 481b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (jc + jc)); 482b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 483b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 484b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org case 4: /* transform for factor of 4 */ 485b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ispan = kspan; 486b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan /= 4; 487b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 488b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 489b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = 1.0; 490b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = 0.0; 491b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 492b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 493b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk + kspan; 494b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = k1 + kspan; 495b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k3 = k2 + kspan; 496b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akp = Re [kk] + Re [k2]; 497b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akm = Re [kk] - Re [k2]; 498b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ajp = Re [k1] + Re [k3]; 499b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ajm = Re [k1] - Re [k3]; 500b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkp = Im [kk] + Im [k2]; 501b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkm = Im [kk] - Im [k2]; 502b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bjp = Im [k1] + Im [k3]; 503b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bjm = Im [k1] - Im [k3]; 504b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] = akp + ajp; 505b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] = bkp + bjp; 506b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ajp = akp - ajp; 507b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bjp = bkp - bjp; 508b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (iSign < 0) { 509b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akp = akm + bjm; 510b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkp = bkm - ajm; 511b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akm -= bjm; 512b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkm += ajm; 513b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } else { 514b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akp = akm - bjm; 515b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkp = bkm + ajm; 516b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akm += bjm; 517b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkm -= ajm; 518b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 519b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* avoid useless multiplies */ 520b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (s1 == 0.0) { 521b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = akp; 522b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = ajp; 523b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k3] = akm; 524b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = bkp; 525b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = bjp; 526b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k3] = bkm; 527b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } else { 528b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = akp * c1 - bkp * s1; 529b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = ajp * c2 - bjp * s2; 530b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k3] = akm * c3 - bkm * s3; 531b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = akp * s1 + bkp * c1; 532b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = ajp * s2 + bjp * c2; 533b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k3] = akm * s3 + bkm * c3; 534b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 535b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k3 + kspan; 536b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < nt); 537b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 538b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 = c1 - (cd * c1 + sd * s1); 539b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = sd * c1 - cd * s1 + s1; 540b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = 2.0 - (c2 * c2 + s1 * s1); 541b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 *= c1; 542b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 *= c2; 543b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* values of c2, c3, s2, s3 that will get used next time */ 544b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 = c1 * c1 - s1 * s1; 545b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s2 = 2.0 * c1 * s1; 546b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c3 = c2 * c1 - s2 * s1; 547b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s3 = c2 * s1 + s2 * c1; 548b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = kk - nt + jc; 549b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < kspan); 550b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = kk - kspan + inc; 551b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < jc); 552b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kspan == jc) 553b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org goto Permute_Results_Label; /* exit infinite loop */ 554b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 555b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 556b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org default: 557b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* transform for odd factors */ 558b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#ifdef FFT_RADIX4 559b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return -1; 560b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 561b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#else /* FFT_RADIX4 */ 562b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = fftstate->factor [ii - 1]; 563b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ispan = kspan; 564b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan /= k; 565b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 566b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org switch (k) { 567b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org case 3: /* transform for factor of 3 (optional code) */ 568b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 569b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 570b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk + kspan; 571b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = k1 + kspan; 572b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = Re [kk]; 573b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = Im [kk]; 574b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aj = Re [k1] + Re [k2]; 575b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bj = Im [k1] + Im [k2]; 576b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] = ak + aj; 577b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] = bk + bj; 578b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak -= 0.5 * aj; 579b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk -= 0.5 * bj; 580b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aj = (Re [k1] - Re [k2]) * s60; 581b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bj = (Im [k1] - Im [k2]) * s60; 582b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = ak - bj; 583b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = ak + bj; 584b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = bk + aj; 585b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = bk - aj; 586b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k2 + kspan; 587b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (nn - 1)); 588b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk -= nn; 589b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < kspan); 590b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 591b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 592b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org case 5: /* transform for factor of 5 (optional code) */ 593b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 = c72 * c72 - s72 * s72; 594b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s2 = 2.0 * c72 * s72; 595b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 596b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 597b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk + kspan; 598b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = k1 + kspan; 599b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k3 = k2 + kspan; 600b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k4 = k3 + kspan; 601b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akp = Re [k1] + Re [k4]; 602b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org akm = Re [k1] - Re [k4]; 603b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkp = Im [k1] + Im [k4]; 604b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bkm = Im [k1] - Im [k4]; 605b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ajp = Re [k2] + Re [k3]; 606b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ajm = Re [k2] - Re [k3]; 607b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bjp = Im [k2] + Im [k3]; 608b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bjm = Im [k2] - Im [k3]; 609b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aa = Re [kk]; 610b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bb = Im [kk]; 611b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] = aa + akp + ajp; 612b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] = bb + bkp + bjp; 613b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = akp * c72 + ajp * c2 + aa; 614b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = bkp * c72 + bjp * c2 + bb; 615b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aj = akm * s72 + ajm * s2; 616b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bj = bkm * s72 + bjm * s2; 617b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = ak - bj; 618b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k4] = ak + bj; 619b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = bk + aj; 620b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k4] = bk - aj; 621b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = akp * c2 + ajp * c72 + aa; 622b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = bkp * c2 + bjp * c72 + bb; 623b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aj = akm * s2 - ajm * s72; 624b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bj = bkm * s2 - bjm * s72; 625b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = ak - bj; 626b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k3] = ak + bj; 627b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = bk + aj; 628b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k3] = bk - aj; 629b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k4 + kspan; 630b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (nn-1)); 631b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk -= nn; 632b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < kspan); 633b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 634b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 635b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org default: 636b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (k != jf) { 637b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jf = k; 638b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = pi2 / (double) k; 639b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = cos(s1); 640b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = sin(s1); 641b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (jf > max_factors){ 642b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return -1; 643b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 644b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Cos [jf - 1] = 1.0; 645b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Sin [jf - 1] = 0.0; 646b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 647b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 648b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1; 649b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1; 650b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k--; 651b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Cos [k - 1] = Cos [j - 1]; 652b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Sin [k - 1] = -Sin [j - 1]; 653b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 654b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j < k); 655b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 656b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 657b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 658b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk; 659b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kk + ispan; 660b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = aa = Re [kk]; 661b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = bb = Im [kk]; 662b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 663b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 += kspan; 664b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 665b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 -= kspan; 666b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 667b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Rtmp [j - 1] = Re [k1] + Re [k2]; 668b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak += Rtmp [j - 1]; 669b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Itmp [j - 1] = Im [k1] + Im [k2]; 670b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk += Itmp [j - 1]; 671b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 672b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Rtmp [j - 1] = Re [k1] - Re [k2]; 673b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Itmp [j - 1] = Im [k1] - Im [k2]; 674b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 += kspan; 675b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k1 < k2); 676b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] = ak; 677b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] = bk; 678b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk; 679b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kk + ispan; 680b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 681b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 682b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 += kspan; 683b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 -= kspan; 684b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj = j; 685b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = aa; 686b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = bb; 687b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aj = 0.0; 688b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bj = 0.0; 689b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = 1; 690b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 691b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k++; 692b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak += Rtmp [k - 1] * Cos [jj - 1]; 693b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk += Itmp [k - 1] * Cos [jj - 1]; 694b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k++; 695b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org aj += Rtmp [k - 1] * Sin [jj - 1]; 696b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bj += Itmp [k - 1] * Sin [jj - 1]; 697b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj += j; 698b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (jj > jf) { 699b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj -= jf; 700b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 701b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k < jf); 702b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = jf - j; 703b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = ak - bj; 704b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = bk + aj; 705b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k2] = ak + bj; 706b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k2] = bk - aj; 707b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 708b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j < k); 709b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += ispan; 710b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < nn); 711b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk -= nn; 712b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < kspan); 713b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 714b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 715b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 716b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* multiply by rotation factor (except for factors of 2 and 4) */ 717b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (ii == mfactor) 718b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org goto Permute_Results_Label; /* exit infinite loop */ 719b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = jc; 720b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 721b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 = 1.0 - cd; 722b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 = sd; 723b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 724b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = c2; 725b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s2 = s1; 726b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += kspan; 727b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 728b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 729b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = Re [kk]; 730b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [kk] = c2 * ak - s2 * Im [kk]; 731b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [kk] = s2 * ak + c2 * Im [kk]; 732b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += ispan; 733b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < nt); 734b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = s1 * s2; 735b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s2 = s1 * c2 + c1 * s2; 736b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 = c1 * c2 - ak; 737b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = kk - nt + kspan; 738b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < ispan); 739b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 = c1 - (cd * c1 + sd * s1); 740b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 += sd * c1 - cd * s1; 741b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c1 = 2.0 - (c2 * c2 + s1 * s1); 742b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org s1 *= c1; 743b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org c2 *= c1; 744b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = kk - ispan + jc; 745b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < kspan); 746b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = kk - kspan + jc + inc; 747b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (jc + jc)); 748b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; 749b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org#endif /* FFT_RADIX4 */ 750b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 751b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 752b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 753b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* permute the results to normal order---done in two stages */ 754b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* permutation for square factors of n */ 755b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.orgPermute_Results_Label: 756b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [0] = ns; 757b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kt) { 758b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = kt + kt + 1; 759b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (mfactor < k) 760b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k--; 761b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 762b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [k] = jc; 763b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 764b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1]; 765b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1]; 766b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 767b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k--; 768b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j < k); 769b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k3 = fftstate->Perm [k]; 770b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan = fftstate->Perm [1]; 771b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = jc; 772b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kspan; 773b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 774b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (nPass != nTotal) { 775b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* permutation for multivariate transform */ 776b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Permute_Multi_Label: 777b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 778b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 779b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = kk + jc; 780b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 781b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ 782b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; 783b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; 784b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += inc; 785b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 += inc; 786b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (k-1)); 787b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += ns - jc; 788b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 += ns - jc; 789b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (nt-1)); 790b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = k2 - nt + kspan; 791b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = kk - nt + jc; 792b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k2 < (ns-1)); 793b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 794b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 795b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 -= fftstate->Perm [j - 1]; 796b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 797b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = fftstate->Perm [j] + k2; 798b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k2 > fftstate->Perm [j - 1]); 799b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 800b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 801b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kk < (k2-1)) 802b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org goto Permute_Multi_Label; 803b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += jc; 804b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 += kspan; 805b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k2 < (ns-1)); 806b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (ns-1)); 807b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } else { 808b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* permutation for single-variate transform (optional code) */ 809b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Permute_Single_Label: 810b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 811b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ 812b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; 813b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; 814b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += inc; 815b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 += kspan; 816b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k2 < (ns-1)); 817b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 818b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 819b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 -= fftstate->Perm [j - 1]; 820b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 821b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = fftstate->Perm [j] + k2; 822b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k2 >= fftstate->Perm [j - 1]); 823b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 1; 824b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 825b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kk < k2) 826b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org goto Permute_Single_Label; 827b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk += inc; 828b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 += kspan; 829b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k2 < (ns-1)); 830b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < (ns-1)); 831b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 832b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jc = k3; 833b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 834b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 835b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if ((kt << 1) + 1 >= mfactor) 836b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return 0; 837b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ispan = fftstate->Perm [kt]; 838b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* permutation for square-free factors of n */ 839b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = mfactor - kt; 840b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [j] = 1; 841b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 842b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->factor [j - 1] *= fftstate->factor [j]; 843b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j--; 844b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j != kt); 845b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kt++; 846b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nn = fftstate->factor [kt - 1] - 1; 847b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (nn > (int) max_perm) { 848b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return -1; 849b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 850b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = jj = 0; 851b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (;;) { 852b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = kt + 1; 853b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = fftstate->factor [kt - 1]; 854b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = fftstate->factor [k - 1]; 855b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 856b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (j > nn) 857b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; /* exit infinite loop */ 858b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj += kk; 859b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org while (jj >= k2) { 860b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj -= k2; 861b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = kk; 862b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k++; 863b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = fftstate->factor [k - 1]; 864b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj += kk; 865b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 866b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [j - 1] = jj; 867b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 868b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* determine the permutation cycles of length greater than 1 */ 869b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = 0; 870b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (;;) { 871b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 872b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j++; 873b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = fftstate->Perm [j - 1]; 874b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk < 0); 875b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (kk != j) { 876b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 877b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = kk; 878b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = fftstate->Perm [k - 1]; 879b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [k - 1] = -kk; 880b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (kk != j); 881b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k3 = kk; 882b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } else { 883b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org fftstate->Perm [j - 1] = -j; 884b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (j == nn) 885b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; /* exit infinite loop */ 886b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 887b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 888b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org max_factors *= inc; 889b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org /* reorder a and b, following the permutation cycles */ 890b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org for (;;) { 891b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j = k3 + 1; 892b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org nt -= ispan; 893b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org ii = nt - inc + 1; 894b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (nt < 0) 895b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org break; /* exit infinite loop */ 896b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 897b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 898b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org j--; 899b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (fftstate->Perm [j - 1] < 0); 900b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj = jc; 901b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 902b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan = jj; 903b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org if (jj > max_factors) { 904b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kspan = max_factors; 905b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 906b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org jj -= kspan; 907b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = fftstate->Perm [j - 1]; 908b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = jc * k + ii + jj; 909b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk + kspan - 1; 910b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = 0; 911b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 912b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2++; 913b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Rtmp [k2 - 1] = Re [k1]; 914b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Itmp [k2 - 1] = Im [k1]; 915b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 -= inc; 916b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k1 != (kk-1)); 917b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 918b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk + kspan - 1; 919b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = k1 - jc * (k + fftstate->Perm [k - 1]); 920b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k = -fftstate->Perm [k - 1]; 921b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 922b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = Re [k2]; 923b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = Im [k2]; 924b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 -= inc; 925b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 -= inc; 926b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k1 != (kk-1)); 927b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org kk = k2 + 1; 928b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k != j); 929b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 = kk + kspan - 1; 930b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2 = 0; 931b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org do { 932b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k2++; 933b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Re [k1] = Rtmp [k2 - 1]; 934b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org Im [k1] = Itmp [k2 - 1]; 935b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org k1 -= inc; 936b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (k1 != (kk-1)); 937b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (jj); 938b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } while (j != 1); 939b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org } 940b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org return 0; /* exit point here */ 941b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org} 942b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org/* ---------------------- end-of-file (c source) ---------------------- */ 943b015cbede88899f67a53fbbe581b02ce8e32794andrew@webrtc.org 944