General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package
Description:
    A package to calculate Discrete Fourier/Cosine/Sine Transforms of 
    1-dimensional sequences of length 2^N.
Files:
    fft4g.c    : FFT Package in C       - Fast Version   I   (radix 4,2)
    fft4g.f    : FFT Package in Fortran - Fast Version   I   (radix 4,2)
    fft4g_h.c  : FFT Package in C       - Simple Version I   (radix 4,2)
    fft8g.c    : FFT Package in C       - Fast Version   II  (radix 8,4,2)
    fft8g.f    : FFT Package in Fortran - Fast Version   II  (radix 8,4,2)
    fft8g_h.c  : FFT Package in C       - Simple Version II  (radix 8,4,2)
    fftsg.c    : FFT Package in C       - Fast Version   III (Split-Radix)
    fftsg.f    : FFT Package in Fortran - Fast Version   III (Split-Radix)
    fftsg_h.c  : FFT Package in C       - Simple Version III (Split-Radix)
    readme.txt : Readme File
    sample1/   : Test Directory
        Makefile    : for gcc, cc
        Makefile.f77: for Fortran
        testxg.c    : Test Program for "fft*g.c"
        testxg.f    : Test Program for "fft*g.f"
        testxg_h.c  : Test Program for "fft*g_h.c"
    sample2/   : Benchmark Directory
        Makefile    : for gcc, cc
        Makefile.pth: POSIX Thread version
        pi_fft.c    : PI(= 3.1415926535897932384626...) Calculation Program
                      for a Benchmark Test for "fft*g.c"
Difference of the Files:
    C and Fortran versions are equal and 
    the same routines are in each version.
    "fft4g*.*" are optimized for most machines.
    "fft8g*.*" are fast on the UltraSPARC.
    "fftsg*.*" are optimized for the machines that 
    have the multi-level (L1,L2,etc) cache.
    The simple versions "fft*g_h.c" use no work area, but 
    the fast versions "fft*g.*" use work areas.
    The fast versions "fft*g.*" have the same specification.
Routines in the Package:
    cdft: Complex Discrete Fourier Transform
    rdft: Real Discrete Fourier Transform
    ddct: Discrete Cosine Transform
    ddst: Discrete Sine Transform
    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
Usage:
    Please refer to the comments in the "fft**.*" file which 
    you want to use. Brief explanations are in the block 
    comments of each package. The examples are also given in 
    the test programs.
Method:
    -------- cdft --------
    fft4g*.*, fft8g*.*:
        A method of in-place, radix 2^M, Sande-Tukey (decimation in 
        frequency). Index of the butterfly loop is in bit 
        reverse order to keep continuous memory access.
    fftsg*.*:
        A method of in-place, Split-Radix, recursive fast 
        algorithm.
    -------- rdft --------
    A method with a following butterfly operation appended to "cdft".
    In forward transform :
        A[k] = sum_j=0^n-1 a[j]*W(n)^(j*k), 0<=k<=n/2, 
            W(n) = exp(2*pi*i/n), 
    this routine makes an array x[] :
        x[j] = a[2*j] + i*a[2*j+1], 0<=j<n/2
    and calls "cdft" of length n/2 :
        X[k] = sum_j=0^n/2-1 x[j] * W(n/2)^(j*k), 0<=k<n.
    The result A[k] are :
        A[k]     = X[k]     - (1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k])), 
        A[n/2-k] = X[n/2-k] + 
                        conjg((1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k]))), 
            0<=k<=n/2
        (notes: conjg() is a complex conjugate, X[n/2]=X[0]).
    -------- ddct --------
    A method with a following butterfly operation appended to "rdft".
    In backward transform :
        C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n, 
    this routine makes an array r[] :
        r[0] = a[0], 
        r[j]   = Re((a[j] - i*a[n-j]) * W(4*n)^j*(1+i)/2), 
        r[n-j] = Im((a[j] - i*a[n-j]) * W(4*n)^j*(1+i)/2), 
            0<j<=n/2
    and calls "rdft" of length n :
        A[k] = sum_j=0^n-1 r[j]*W(n)^(j*k), 0<=k<=n/2, 
            W(n) = exp(2*pi*i/n).
    The result C[k] are :
        C[2*k]   =  Re(A[k] * (1-i)), 
        C[2*k-1] = -Im(A[k] * (1-i)).
    -------- ddst --------
    A method with a following butterfly operation appended to "rdft".
    In backward transform :
        S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n, 
    this routine makes an array r[] :
        r[0] = a[0], 
        r[j]   = Im((a[n-j] - i*a[j]) * W(4*n)^j*(1+i)/2), 
        r[n-j] = Re((a[n-j] - i*a[j]) * W(4*n)^j*(1+i)/2), 
            0<j<=n/2
    and calls "rdft" of length n :
        A[k] = sum_j=0^n-1 r[j]*W(n)^(j*k), 0<=k<=n/2, 
            W(n) = exp(2*pi*i/n).
    The result S[k] are :
        S[2*k]   =  Re(A[k] * (1+i)), 
        S[2*k-1] = -Im(A[k] * (1+i)).
    -------- dfct --------
    A method to split into "dfct" and "ddct" of half length.
    The transform :
        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
    is divided into :
        C[2*k]   = sum'_j=0^n/2  (a[j]+a[n-j])*cos(pi*j*k/(n/2)), 
        C[2*k+1] = sum_j=0^n/2-1 (a[j]-a[n-j])*cos(pi*j*(k+1/2)/(n/2))
        (sum' is a summation whose last term multiplies 1/2).
    This routine uses "ddct" recursively.
    To keep the in-place operation, the data in fft*g_h.*
    are sorted in bit reversal order.
    -------- dfst --------
    A method to split into "dfst" and "ddst" of half length.
    The transform :
        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
    is divided into :
        S[2*k]   = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(pi*j*k/(n/2)), 
        S[2*k+1] = sum'_j=1^n/2  (a[j]+a[n-j])*sin(pi*j*(k+1/2)/(n/2))
        (sum' is a summation whose last term multiplies 1/2).
    This routine uses "ddst" recursively.
    To keep the in-place operation, the data in fft*g_h.*
    are sorted in bit reversal order.
Reference:
    * Masatake MORI, Makoto NATORI, Tatuo TORII: Suchikeisan, 
      Iwanamikouzajyouhoukagaku18, Iwanami, 1982 (Japanese)
    * Henri J. Nussbaumer: Fast Fourier Transform and Convolution 
      Algorithms, Springer Verlag, 1982
    * C. S. Burrus, Notes on the FFT (with large FFT paper list)
      http://www-dsp.rice.edu/research/fft/fftnote.asc
    email: ooura@mmm.t.u-tokyo.ac.jp
    download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
    You may use, copy, modify this code for any purpose and 
    without fee. You may distribute this ORIGINAL package.
History:
    Dec. 1995  : Edit the General Purpose FFT
    Mar. 1996  : Change the specification
    Jun. 1996  : Change the method of trigonometric function table
    Sep. 1996  : Modify the documents
    Feb. 1997  : Change the butterfly loops
    Dec. 1997  : Modify the documents
    Dec. 1997  : Add "fft4g.*"
    Jul. 1998  : Fix some bugs in the documents
    Jul. 1998  : Add "fft8g.*" and delete "fft4f.*"
    Jul. 1998  : Add a benchmark program "pi_fft.c"
    Jul. 1999  : Add a simple version "fft*g_h.c"
    Jul. 1999  : Add a Split-Radix FFT package "fftsg*.c"
    Sep. 1999  : Reduce the memory operation (minor optimization)
    Oct. 1999  : Change the butterfly structure of "fftsg*.c"
    Oct. 1999  : Save the code size
    Sep. 2001  : Add "fftsg.f"
    Sep. 2001  : Add Pthread & Win32thread routines to "fftsg*.c"
    Dec. 2006  : Fix a minor bug in "fftsg.f"