* * * * * * * * * * * * * * * * * VFFTPK * * * * * * * * * * * * * * * * * (Version 2.1 May 1990) A vectorized package of Fortran subprograms for the fast Fourier transform of multiple real sequences by Roland A. Sweet and Linda L. Lindgren Center for Computing and Applied Mathematics National Institute of Standards and Technology Boulder, Colorado 80303 Ronald F. Boisvert Center for Computing and Applied Mathematics National Institute of Standards and Technology Gaithersburg, Maryland 20899 * * * * * * * * * * * * * * * * * Background * * * * * * * * * * * * * * * * * NOTATION -------- To simulate the mathematical symbol 'sigma' that represents summation, we define the notation SUM(l=is,if)[ y(l) ] = y(is) + y(is+1) + . . . + y(if) . DISCRETE FOURIER TRANSFORM -------------------------- Given a complex periodic sequence x(j), j=0,1,...,n-1, the discrete Fourier transform is defined by x(j) = sqrt(1/n)*SUM(k=0,n-1)[ f(k)*exp(2j*k*i*pi/n) ] , for j = 0, 1, . . . , n-1 , where i = sqrt(-1). The discrete Fourier coefficients, f(k), are given by f(k) = sqrt(1/n)*SUM(j=0,n-1)[ x(j)*exp(-2j*k*i*pi/n) ] , for k = 0, 1, . . . , n-1 . REAL PERIODIC TRANSFORM ----------------------- To facilitate the following discussion we define the integer m by m = n/2 when n is even, and m = (n+1)/2 when n is odd. If the given sequence x(j) is real, the Fourier coefficients are conjugate symmetric, i.e. f(n-k) is equal to the complex conjugate of f(k). This relation implies that the sequence x(j) can be completely determined from the n real quantities real(f(0)) = sqrt(1/n)*SUM(j=0,n-1)[ x(j) ], real(f(k)) = sqrt(1/n)*SUM(j=0,n-1)[ x(j)*cos(2j*k*pi/n) ] imag(f(k)) = -sqrt(1/n)*SUM(j=0,n-1)[ x(j)*sin(2j*k*pi/n) ] for k = 1, 2, . . . , m-1, and, when n is even, real(f(n/2)) = sqrt(1/n)*SUM(j=0,n-1)[ (-1)**j*x(j) ]. TRIGONOMETRIC REPRESENTATION ---------------------------- When the sequence x(j) is real the discrete Fourier transform may be re-written as the trigonometric series x(j) = sqrt(2/n)* c(0)/2 + (-1)**j*c(n-1)/2 + SUM(k=1,m-1) [c(2k-1)*cos(2j*k*pi/n) + c(2k)*sin(2j*k*pi/n)] when n is an even integer (n = 2*m), or as x(j) = sqrt(2/n)* c(0)/2 + SUM(k=1,m-1) [c(2k-1)*cos(2j*k*pi/n) + c(2k)*sin(2j*k*pi/n)] when n is an odd integer (n = 2*m-1), for j = 0, 1, . . . , n-1 . The coefficients c(k) are defined by the relations c(0) = sqrt(2/n)*SUM(j=0,n-1)[ x(j) ] for k = 1, 2, . . . , m-1 , c(2*k-1) = sqrt(2/n)*SUM(j=0,n-1)[ x(j)*cos(2j*k*pi/n) ] c(2*k) = sqrt(2/n)*SUM(j=0,n-1)[ x(j)*sin(2j*k*pi/n) ] and, when n is even c(n-1) = sqrt(2/n)*SUM(j=0,n-1)[ (-1)**j*x(j) ] . The relationship between the coefficients c(k) and the conjugate symmetric Fourier coefficients f(k) is c(0) = sqrt(2)*f(0) , c(2*k-1) = sqrt(2)*real(f(k)) and c(2*k) = -sqrt(2)*imag(f(k)) for k = 1, 2, . . . , m-1, and when n is even c(n-1) = sqrt(2)*f(n/2) . SINE TRANSFORM OF AN ODD SEQUENCE --------------------------------- Let x(k) be a real odd sequence, that is, it can be expanded in terms of a trigonometric series that contains only sine terms. The discrete odd Fourier transform is given by c(k) = SUM(j=1,n-1)[ 2x(j)*sin(j*k*pi/n) ]/sqrt(2n). for k = 1, 2, . . ., n-1. The inverse transform is identical, that is, x(k) = SUM(j=1,n-1)[ 2c(j)*sin(j*k*pi/n) ]/sqrt(2n) for k = 1, 2, . . ., n-1. COSINE TRANSFORM OF AN EVEN SEQUENCE ------------------------------------ Let x(k) be a real even sequence, that is, it can be expanded in terms of a trigonometric series that contains only cosine terms. The discrete even Fourier transform is given by c(k) = ( x(0) + (-1)**k*x(n) + SUM(j=1,n-1)[ 2x(j)*cos(j*k*pi/n) ] )/sqrt(2n) for k = 0, 1, 2, . . ., n. The inverse transform is identical, that is, x(k) = ( c(0) + (-1)**k*c(n) + SUM(j=1,n-1)[ 2c(j)*cos(j*k*pi/n) ] )/sqrt(2n) for k = 0, 1, 2, . . ., n. QUARTER-WAVE COSINE TRANSFORM ----------------------------- Let x(k) be a real even quarter-wave sequence, that is, it can be expanded in terms of a cosine series with only odd wave numbers. The discrete cosine quarter-wave Fourier transform is given by c(k) = ( x(1) + SUM(j=2,n)[ 2x(j)*cos((j-1)*(2k-1)*pi/2n) ] )/sqrt(4n) for k = 1, 2, . . ., n. The inverse transform is given by x(k) = SUM(j=1,n)[ 4c(j)*cos((2j-1)*(k-1)*pi/2n) ]/sqrt(4n) for k = 1, 2, . . ., n. QUARTER-WAVE SINE TRANSFORM --------------------------- Let x(k) be a real odd quarter-wave sequence; that is, it can be expanded in terms of a sine series with only odd wave numbers. The discrete quarter-wave sine transform is given by c(k) = ( SUM(j=1,n-1)[ 2x(j)*sin(j*(2k-1)*pi/2n) ] + (-1)**(k-1)*x(n) )/sqrt(4n) for k = 1, 2, . . ., n. The inverse transform is given by x(k) = ( SUM(j=1,n) [ 4c(j)*sin((2j-1)*k*pi/2n) ] )/sqrt(4n) for k = 1, 2, . . ., n. * * * * * * * * * * * * * * * * * * * * * * * * Description of Package * * * * * * * * * * * * * * * * * * * * * * * * This package consists of subroutines that compute the transforms between sequences x(j) and the Fourier coefficients f(k) described above for multiple real sequences. The real periodic transform is a variant of the Stockham autosort algorithm. There is no restriction on the length of the sequence. The other specialized transforms are implemented by preprocessing and postprocessing. In the description below we use the following terms: forward transform = Fourier analysis. Computation of Fourier coefficients f(k) from the sequence x(k). backward transform = Fourier synthesis. Computation of the sequence x(k) from its Fourier coefficients f(k). The user-callable subroutines are : 1. VRFFTI Initialization routine for VRFFTF and VRFFTB 2. VRFFTF Forward transform of multiple real periodic sequences 3. VRFFTB Backward transform for multiple real periodic sequences 4. VSINTI Initialization routine for VSINT 5. VSINT Forward or backward transform for multiple odd sequences 6. VCOSTI Initialization routine for VCOST 7. VCOST Forward or backward transform for multiple even sequences 8. VSINQI Initialization routine for VSINQF and VSINQB 9. VSINQF Forward transform of multiple odd sequences with only odd wave numbers (quarter-wave sine transform) 10. VSINQB Backward transform for multiple odd sequences with only odd wave numbers (quarter-wave sine transform) 11 VCOSQI Initialization routine for VCOSQF and VCOSQB 12 VCOSQF Forward transform of multiple even sequences with only odd wave numbers (quarter-wave cosine transform) 13. VCOSQB Backward transform of multiple even sequences with only odd wave numbers (quarter-wave cosine transform) This package is a straightforward vectorization of the scalar package FFTPACK (Version 3, June 1979) written by Paul N. Swarztrauber of the National Center for Atmospheric Research, Boulder, Colorado, 80307. All vector lengths are equal to the number of sequences being transformed. One slight change has been made: symmetric scaling has been incorporated into the forward and backward transforms. Revision history: Version 1.0 (Aug 1985) : VRFFTF, VRFFTB, VRFFTI developed by Sweet and Lindgren Version 2.0 (Mar 1987) : VSINT, VSINTI, VCOST, VCOSTI, VSINQF, VSINQB, VSINQI, VCOSQF, VCOSQB, VCOSQI added by Boisvert Version 2.1 (May 1990) : documentation revised To get a particular transform, use a request like send vrffti vrfftf from vfftpk. To get the entire library, say send all from vfftpk.