 
 
                     * * * * * * * * * * * * * *
                     *                         *
                     *          VFFTPK         *
                     *                         *
                     * * * * * * * * * * * * * *
 
                       (Version 2  March 1987)
 
          A vectorized package of Fortran subprograms for the
      fast Fourier transform of multiple real periodic sequences
 
                                  by
 
                Roland A. Sweet and Linda L. Lindgren
 
                    Scientific Computing Division
                     National Bureau of Standards
                       Boulder, Colorado 80303
 
                         Ronald F. Boisvert
 
                    Scientific Computing Division
                     National Bureau of Standards
                     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 vectors being
 transformed. One slight change has been made; symmetric scaling has
 been incorporated into the forward and backward transforms.
 
 
