fftxf.doc

FFT代原碼為C++需要測過才能用

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]

    fft2f.c    : FFT Package in C       - Easy Version
    fft2f.f    : FFT Package in Fortran - Easy Version
    fft2ft.c   : Test Program of "fft2f.c"
    fft2ft.f   : Test Program of "fft2f.f"
    fft4f.c    : FFT Package in C       - Fast Version
    fft4f.f    : FFT Package in Fortran - Fast Version
    fft4ft.c   : Test Program of "fft4f.c"
    fft4ft.f   : Test Program of "fft4f.f"
    fftxf.doc  : Document of "fft*f.*"
    fftxfj.doc : Document of "fft*f.*" (Japanese)


[Difference of Files]

    C and Fortran versions are same.
    Easy versions use no work area.
    Fast versions use work areas.


[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]

    Comments :
        Brief explanations are in block comments of each packages.
        The examples are given in the test programs.
        The cos/sin table and work areas are not necessary if you use 
        "fft2f.*".

    -------- Complex DFT (Discrete Fourier Transform) --------
        <CASE1>
            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
        <CASE2>
            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)

        [usage]
            in fft2f.c
                <CASE1>
                    cdft(2*n, cos(M_PI/n), sin(M_PI/n), a);
                <CASE2>
                    cdft(2*n, cos(M_PI/n), -sin(M_PI/n), a);
            in fft2f.f
                <CASE1>
                    call cdft(2*n, cos(pi/n), sin(pi/n), a)
                <CASE2>
                    call cdft(2*n, cos(pi/n), -sin(pi/n), a)
            in fft4f.c
                <CASE1>
                    ip[0] = 0; // first time only
                    cdft(2*n, 1, a, ip, w);
                <CASE2>
                    ip[0] = 0; // first time only
                    cdft(2*n, -1, a, ip, w);
            in fft4f.f
                <CASE1>
                    ip(0) = 0  ! first time only
                    call cdft(2*n, 1, a, ip, w)
                <CASE2>
                    ip(0) = 0  ! first time only
                    call cdft(2*n, -1, a, ip, w)

        [arguments]
            2*n : real data length (n is a complex data length)
                      type   : integer
                      detail : n >= 1, n = power of 2
            a   : input/output data
                      type   : double real array
                      size   : from 0 to 2*n-1
                      detail : input data
                                   a[2*j] = Re(x[j]), 
                                   a[2*j+1] = Im(x[j]), 0<=j<n
                               output data
                                   a[2*k] = Re(X[k]), 
                                   a[2*k+1] = Im(X[k]), 0<=k<n

            ip  : work area for bit reversal
                      type   : integer array
                      size   : from 0 to 1+sqrt(n)
                               strictly, 
                               length of ip >= 2+sqrt(n)  ; if n mod 4 == 0
                                               2+sqrt(n/2); otherwise
                      detail : ip[2...*] is the work area, ip[0],ip[1] are 
                               pointers of the cos/sin table.
            w   : cos/sin table
                      type   : double real array
                      size   : from 0 to n/2-1
                      detail : w[] is calculated if ip[0] == 0.
                               If you call this function repeatedly, the 
                               recalculation of the table is not necessary.

        [caution!]
            Conventional definitions of DFT / Inverse of DFT are :
                a) CASE1 defines backward, 1/n scaled CASE2 defines forward
                b) CASE1 defines forward, 1/n scaled CASE2 defines backward
                c) 1/n scaled CASE1 defines backward, CASE2 defines forward
                d) 1/n scaled CASE1 defines forward, CASE2 defines backward
                e) etc.
            You must know which one you use.

    -------- Real DFT / Inverse of Real DFT --------
        <CASE1> (RDFT)
            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
        <CASE2> (IRDFT)
            A[k] = R[0]/2 + R[n/2]/2 + 
                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n

        [usage]
            in fft2f.c
                <CASE1>
                    rdft(n, cos(M_PI/n), sin(M_PI/n), a);
                <CASE2>
                    rdft(n, cos(M_PI/n), -sin(M_PI/n), a);
            in fft2f.f
                <CASE1>
                    call rdft(n, cos(pi/n), sin(pi/n), a)
                <CASE2>
                    call rdft(n, cos(pi/n), -sin(pi/n), a)
            in fft4f.c
                <CASE1>
                    ip[0] = 0; // first time only
                    rdft(n, 1, a, ip, w);
                <CASE2>
                    ip[0] = 0; // first time only
                    rdft(n, -1, a, ip, w);
            in fft4f.f
                <CASE1>
                    ip(0) = 0  ! first time only
                    call rdft(n, 1, a, ip, w)
                <CASE2>
                    ip(0) = 0  ! first time only
                    call rdft(n, -1, a, ip, w)

        [arguments]
            n   : real data length
                      type   : integer
                      detail : n >= 2, n = power of 2
            a   : input/output data
                      type   : double real array
                      size   : from 0 to n-1
                      detail : output data in CASE1
                                   a[2*k] = R[k], 0<=k<n/2
                                   a[2*k+1] = I[k], 0<k<n/2
                                   a[1] = R[n/2]
                               input data in CASE2
                                   a[2*j] = R[j], 0<=j<n/2
                                   a[2*j+1] = I[j], 0<j<n/2
                                   a[1] = R[n/2]

            ip  : work area for bit reversal
                      type   : integer array
                      size   : from 0 to 1+sqrt(n/2)
                               strictly, 
                               length of ip >= 2+sqrt(n/2); if n mod 4 == 2
                                               2+sqrt(n/4); otherwise
                      detail : ip[2...*] is the work area, ip[0],ip[1] are 
                               pointers of the cos/sin table.
            w   : cos/sin table
                      type   : double real array
                      size   : from 0 to n/2-1
                      detail : w[] is calculated if ip[0] == 0.
                               If you call this function repeatedly, the 
                               recalculation of the table is not necessary.

        [remark]
            Inverse of the transform CASE1 is equal to the transform CASE2 
            with scaling :
                a[j] = a[j] * 2/n, 0<=j<n.

    -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
        <CASE1> (IDCT)
            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
        <CASE2> (DCT)
            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n

        [usage]
            in fft2f.c
                <CASE1>
                    ddct(n, cos(M_PI/n/2), sin(M_PI/n/2), a);
                <CASE2>
                    ddct(n, cos(M_PI/n/2), -sin(M_PI/n/2), a);
            in fft2f.f
                <CASE1>
                    call ddct(n, cos(pi/n/2), sin(pi/n/2), a)
                <CASE2>
                    call ddct(n, cos(pi/n/2), -sin(pi/n/2), a)
            in fft4f.c
                <CASE1>
                    ip[0] = 0; // first time only
                    ddct(n, 1, a, ip, w);
                <CASE2>
                    ip[0] = 0; // first time only
                    ddct(n, -1, a, ip, w);
            in fft4f.f
                <CASE1>
                    ip(0) = 0  ! first time only
                    call ddct(n, 1, a, ip, w)
                <CASE2>
                    ip(0) = 0  ! first time only
                    call ddct(n, -1, a, ip, w)

        [arguments]
            n   : data length
                      type   : integer
                      detail : n >= 2, n = power of 2
            a   : input/output data
                      type   : double real array
                      size   : from 0 to n-1
                      detail : output data
                                   a[k] = C[k], 0<=k<n

            ip  : work area for bit reversal
                      type   : integer array
                      size   : from 0 to 1+sqrt(n/2)
                               strictly, 
                               length of ip >= 2+sqrt(n/2); if n mod 4 == 2
                                               2+sqrt(n/4); otherwise
                      detail : ip[2...*] is the work area, ip[0],ip[1] are 
                               pointers of the cos/sin table.
            w   : cos/sin table
                      type   : double real array
                      size   : from 0 to n*5/4-1
                      detail : w[] is calculated if ip[0] == 0.
                               If you call this function repeatedly, the 
                               recalculation of the table is not necessary.

        [remark]
            Inverse of the transform CASE2 is equal to the transform CASE1 
            with scaling :
                a[0] = a[0] * 1/n, a[j] = a[j] * 2/n, 0<j<n.

    -------- DST (Discrete Sine Transform) / Inverse of DST --------
        <CASE1> (IDST)
            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
        <CASE2> (DST)
            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n

        [usage]
            in fft2f.c
                <CASE1>
                    ddst(n, cos(M_PI/n/2), sin(M_PI/n/2), a);
                <CASE2>
                    ddst(n, cos(M_PI/n/2), -sin(M_PI/n/2), a);
            in fft2f.f
                <CASE1>
                    call ddst(n, cos(pi/n/2), sin(pi/n/2), a)
                <CASE2>
                    call ddst(n, cos(pi/n/2), -sin(pi/n/2), a)
            in fft4f.c
                <CASE1>
                    ip[0] = 0; // first time only
                    ddst(n, 1, a, ip, w);
                <CASE2>
                    ip[0] = 0; // first time only