fft.doc

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                  version 3  june 1979
  a package of fortran subprograms for the fast fourier
      transform of periodic and other symmetric sequences
  paul n swarztrauber
  national center for atmospheric research  boulder,colorado 80307
   which is sponsored by the national science foundation
  modified by P. Bjorstad

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

1.   drffti    initialize  drfftf and drfftb
2.   drfftf    forward transform of a real periodic sequence
3.   drfftb    backward transform of a real coefficient array

4.   defftf    a simplified real periodic forward transform
5.   defftb    a simplified real periodic backward transform

6.   dsinti    initialize dsint
7.   dsint     sine transform of a real odd sequence

8.   dcosti    initialize dcost
9.   dcost     cosine transform of a real even sequence

10.  dsinqi    initialize dsinqf and dsinqb
11.  dsinqf    forward sine transform with odd wave numbers
12.  dsinqb    unnormalized inverse of dsinqf

13.  dcosqi    initialize dcosqf and dcosqb
14.  dcosqf    forward cosine transform with odd wave numbers
15.  dcosqb    unnormalized inverse of dcosqf

16.  dcffti     initialize dcfftf and dcfftb
17.  dcfftf     forward transform of a complex periodic sequence
18.  dcfftb     unnormalized inverse of dcfftf

For single precision, change initial d to s.
Each subroutine is described below. the names used refer to
the double precision version, but the same description
applies for the single precision version.

****************************************************************

subroutine drffti(n,wsave)

****************************************************************

subroutine drffti initializes the array wsave which is used in
both drfftf and drfftb. the prime factorization of n together with
a tabulation of the trigonometric functions are computed and
stored in wsave.

input parameter

n       the length of the sequence to be transformed.

output parameter

wsave   a work array which must be dimensioned at least 2*n+15.
        the same work array can be used for both drfftf and drfftb
        as long as n remains unchanged. different wsave arrays
        are required for different values of n. the contents of
        wsave must not be changed between calls of drfftf or drfftb.

******************************************************************

subroutine drfftf(n,r,wsave)

******************************************************************

subroutine drfftf computes the fourier coefficients of a real
perodic sequence (fourier analysis). the transform is defined
below at output parameter r.

input parameters

n       the length of the array r to be transformed.  the method
        is most efficient when n is a product of small primes.
        n may change so long as different work arrays are provided

r       a real array of length n which contains the sequence
        to be transformed

wsave   a work array which must be dimensioned at least 2*n+15.
        in the program that calls drfftf. the wsave array must be
        initialized by calling subroutine drffti(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.
        the same wsave array can be used by drfftf and drfftb.


output parameters

r       r(1) = the sum from i=1 to i=n of r(i)

        if n is even set l =n/2   , if n is odd set l = (n+1)/2

          then for k = 2,...,l

             r(2*k-2) = the sum from i = 1 to i = n of

                  r(i)*cos((k-1)*(i-1)*2*pi/n)

             r(2*k-1) = the sum from i = 1 to i = n of

                 -r(i)*sin((k-1)*(i-1)*2*pi/n)

        if n is even

             r(n) = the sum from i = 1 to i = n of

                  (-1)**(i-1)*r(i)

 *****  note
             this transform is unnormalized since a call of drfftf
             followed by a call of drfftb will multiply the input
             sequence by n.

wsave   contains results which must not be destroyed between
        calls of drfftf or drfftb.


******************************************************************

subroutine drfftb(n,r,wsave)

******************************************************************

subroutine drfftb computes the real perodic sequence from its
fourier coefficients (fourier synthesis). the transform is defined
below at output parameter r.

input parameters

n       the length of the array r to be transformed.  the method
        is most efficient when n is a product of small primes.
        n may change so long as different work arrays are provided

r       a real array of length n which contains the sequence
        to be transformed

wsave   a work array which must be dimensioned at least 2*n+15.
        in the program that calls drfftb. the wsave array must be
        initialized by calling subroutine drffti(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.
        the same wsave array can be used by drfftf and drfftb.


output parameters

r       for n even and for i = 1,...,n

             r(i) = r(1)+(-1)**(i-1)*r(n)

                  plus the sum from k=2 to k=n/2 of

                   2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)

                  -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)

        for n odd and for i = 1,...,n

             r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of

                  2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)

                 -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)

 *****  note
             this transform is unnormalized since a call of drfftf
             followed by a call of drfftb will multiply the input
             sequence by n.

wsave   contains results which must not be destroyed between
        calls of drfftb or drfftf.


******************************************************************

subroutine defftf(n,r,azero,a,b,wsave)

******************************************************************

subroutine defftf computes the fourier coefficients of a real
perodic sequence (fourier analysis). the transform is defined
below at output parameters azero,a and b. defftf is a simplified
version of drfftf. it is not as fast as drfftf since scaling
and initialization are computed for each transform. the repeated
initialization can be suppressed by removeing the statment
( call deffti(n,wsave) ) from both defftf and defftb and inserting
it at the appropriate place in your program.

input parameters

n       the length of the array r to be transformed.  the method
        is must efficient when n is the product of small primes.

r       a real array of length n which contains the sequence
        to be transformed. r is not destroyed.

wsave   a work array with at least 3*n+15 locations.

output parameters

azero   the sum from i=1 to i=n of r(i)/n

a,b     for n even b(n/2)=0. and a(n/2) is the sum from i=1 to
        i=n of (-1)**(i-1)*r(i)/n

        for n even define kmax=n/2-1
        for n odd  define kmax=(n-1)/2

        then for  k=1,...,kmax

             a(k) equals the sum from i=1 to i=n of

                  2./n*r(i)*cos(k*(i-1)*2*pi/n)

             b(k) equals the sum from i=1 to i=n of

                  2./n*r(i)*sin(k*(i-1)*2*pi/n)


******************************************************************

subroutine defftb(n,r,azero,a,b,wsave)

******************************************************************

subroutine defftb computes a real perodic sequence from its
fourier coefficients (fourier synthesis). the transform is
defined below at output parameter r. defftb is a simplified
version of drfftb. it is not as fast as drfftb since scaling and
initialization are computed for each transform. the repeated
initialization can be suppressed by removeing the statment
( call deffti(n,wsave) ) from both defftf and defftb and inserting
( call deffti(n,wsave) ) from both defftf and defftb and inserting
it at the appropriate place in your program.

input parameters

n       the length of the output array r.  the method is most
        efficient when n is the product of small primes.

azero   the constant fourier coefficient

a,b     arrays which contain the remaining fourier coefficients
        these arrays are not destroyed.

        the length of these arrays depends on whether n is even or
        odd.

        if n is even n/2    locations are required
        if n is odd (n-1)/2 locations are required

wsave   a work array with at least 3*n+15 locations.


output parameters

r       if n is even define kmax=n/2
        if n is odd  define kmax=(n-1)/2

        then for i=1,...,n

             r(i)=azero plus the sum from k=1 to k=kmax of

             a(k)*cos(k*(i-1)*2*pi/n)+b(k)*sin(k*(i-1)*2*pi/n)

********************* complex notation **************************

        for j=1,...,n

        r(j) equals the sum from k=-kmax to k=kmax of

             c(k)*exp(i*k*(j-1)*2*pi/n)

        where

             c(k) = .5*cmplx(a(k),-b(k))   for k=1,...,kmax

             c(-k) = conjg(c(k))

             c(0) = azero

                  and i=sqrt(-1)

*************** amplitude - phase notation ***********************

        for i=1,...,n

        r(i) equals azero plus the sum from k=1 to k=kmax of

             alpha(k)*cos(k*(i-1)*2*pi/n+beta(k))

        where

             alpha(k) = sqrt(a(k)*a(k)+b(k)*b(k))

             cos(beta(k))=a(k)/alpha(k)

             sin(beta(k))=-b(k)/alpha(k)

******************************************************************

subroutine dsinti(n,wsave)

******************************************************************

subroutine dsinti initializes the array wsave which is used in
subroutine dsint. the prime factorization of n together with
a tabulation of the trigonometric functions are computed and
stored in wsave.

input parameter

n       the length of the sequence to be transformed.  the method
        is most efficient when n+1 is a product of small primes.

output parameter

wsave   a work array with at least int(2.5*n+15) locations.
        different wsave arrays are required for different values
        of n. the contents of wsave must not be changed between
        calls of dsint.

******************************************************************

subroutine dsint(n,x,wsave)

******************************************************************

subroutine dsint computes the discrete fourier sine transform
of an odd sequence x(i). the transform is defined below at
output parameter x.

dsint is the unnormalized inverse of itself since a call of dsint
followed by another call of dsint will multiply the input sequence
x by 2*(n+1).

the array wsave which is used by subroutine dsint must be
initialized by calling subroutine dsinti(n,wsave).

input parameters

n       the length of the sequence to be transformed.  the method
        is most efficient when n+1 is the product of small primes.

x       an array which contains the sequence to be transformed

              ************important*************

              x must be dimensioned at least n+1

wsave   a work array with dimension at least int(2.5*n+15)
        in the program that calls dsint. the wsave array must be
        initialized by calling subroutine dsinti(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.

output parameters

x       for i=1,...,n

             x(i)= the sum from k=1 to k=n

                  2*x(k)*sin(k*i*pi/(n+1))

             a call of dsint followed by another call of
             dsint will multiply the sequence x by 2*(n+1).
             hence dsint is the unnormalized inverse
             of itself.

wsave   contains initialization calculations which must not be
        destroyed between calls of dsint.

******************************************************************

subroutine dcosti(n,wsave)

******************************************************************

subroutine dcosti initializes the array wsave which is used in
subroutine dcost. the prime factorization of n together with
a tabulation of the trigonometric functions are computed and
stored in wsave.

input parameter

n       the length of the sequence to be transformed.  the method
        is most efficient when n-1 is a product of small primes.

output parameter

wsave   a work array which must be dimensioned at least 3*n+15.
        different wsave arrays are required for different values
        of n. the contents of wsave must not be changed between
        calls of dcost.

******************************************************************

subroutine dcost(n,x,wsave)

******************************************************************

subroutine dcost computes the discrete fourier cosine transform
of an even sequence x(i). the transform is defined below at output
parameter x.

dcost is the unnormalized inverse of itself since a call of dcost
followed by another call of dcost will multiply the input sequence
x by 2*(n-1). the transform is defined below at output parameter x

the array wsave which is used by subroutine dcost must be
initialized by calling subroutine dcosti(n,wsave).

input parameters

n       the length of the sequence x. n must be greater than 1.
        the method is most efficient when n-1 is a product of
        small primes.

x       an array which contains the sequence to be transformed

wsave   a work array which must be dimensioned at least 3*n+15
        in the program that calls dcost. the wsave array must be
        initialized by calling subroutine dcosti(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.

output parameters

x       for i=1,...,n

           x(i) = x(1)+(-1)**(i-1)*x(n)

             + the sum from k=2 to k=n-1

               2*x(k)*cos((k-1)*(i-1)*pi/(n-1))

             a call of dcost followed by another call of
             dcost will multiply the sequence x by 2*(n-1)
             hence dcost is the unnormalized inverse
             of itself.

wsave   contains initialization calculations which must not be
        destroyed between calls of dcost.

******************************************************************

subroutine dsinqi(n,wsave)

******************************************************************

subroutine dsinqi initializes the array wsave which is used in
both dsinqf and dsinqb. the prime factorization of n together with
a tabulation of the trigonometric functions are computed and
stored in wsave.

input parameter

n       the length of the sequence to be transformed. the method
        is most efficient when n is a product of small primes.

output parameter

wsave   a work array which must be dimensioned at least 3*n+15.
        the same work array can be used for both dsinqf and dsinqb
        as long as n remains unchanged. different wsave arrays
        are required for different values of n. the contents of
        wsave must not be changed between calls of dsinqf or
        dsinqb.

******************************************************************

subroutine dsinqf(n,x,wsave)

******************************************************************

subroutine dsinqf computes the fast fourier transform of quarter
wave data. that is , dsinqf computes the coefficients in a sine
series representation with only odd wave numbers. the transform
is defined below at output parameter x.

dsinqb is the unnormalized inverse of dsinqf -- a call of dsinqf
followed by a call of dsinqb will multiply the input sequence x
by 4*n.

the array wsave which is used by subroutine dsinqf must be
initialized by calling subroutine dsinqi(n,wsave).


input parameters

n       the length of the array x to be transformed.  the method
        is most efficient when n is a product of small primes.

x       an array which contains the sequence to be transformed

wsave   a work array which must be dimensioned at least 3*n+15.
        in the program that calls dsinqf. the wsave array must be
        initialized by calling subroutine dsinqi(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.

output parameters

x       for i=1,...,n

             x(i) = (-1)**(i-1)*x(n)

                + the sum from k=1 to k=n-1 of

                2*x(k)*sin((2*i-1)*k*pi/(2*n))

             a call of dsinqf followed by a call of
             dsinqb will multiply the sequence x by 4*n.
             therefore dsinqb is the unnormalized inverse
             of dsinqf.

wsave   contains initialization calculations which must not
        be destroyed between calls of dsinqf or dsinqb.

******************************************************************

subroutine dsinqb(n,x,wsave)

******************************************************************

subroutine dsinqb computes the fast fourier transform of quarter
wave data. that is , dsinqb computes a sequence from its
representation in terms of a sine series with odd wave numbers.
the transform is defined below at output parameter x.

dsinqf is the unnormalized inverse of dsinqb -- a call of dsinqb
followed by a call of dsinqf will multiply the input sequence x
by 4*n.

the array wsave which is used by subroutine dsinqb must be
initialized by calling subroutine dsinqi(n,wsave).


input parameters

n       the length of the array x to be transformed.  the method
        is most efficient when n is a product of small primes.

x       an array which contains the sequence to be transformed

wsave   a work array which must be dimensioned at least 3*n+15.
        in the program that calls dsinqb. the wsave array must be
        initialized by calling subroutine dsinqi(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.

output parameters

x       for i=1,...,n

             x(i)= the sum from k=1 to k=n of

               4*x(k)*sin((2k-1)*i*pi/(2*n))

             a call of dsinqb followed by a call of
             dsinqf will multiply the sequence x by 4*n.
             therefore dsinqf is the unnormalized inverse
             of dsinqb.

wsave   contains initialization calculations which must not
        be destroyed between calls of dsinqb or dsinqf.

******************************************************************

subroutine dcosqi(n,wsave)

******************************************************************

subroutine dcosqi initializes the array wsave which is used in
both dcosqf and dcosqb. the prime factorization of n together with
a tabulation of the trigonometric functions are computed and
stored in wsave.

input parameter

n       the length of the array to be transformed.  the method
        is most efficient when n is a product of small primes.

output parameter

wsave   a work array which must be dimensioned at least 3*n+15.
        the same work array can be used for both dcosqf and dcosqb
        as long as n remains unchanged. different wsave arrays
        are required for different values of n. the contents of
        wsave must not be changed between calls of dcosqf or
        dcosqb.

******************************************************************

subroutine dcosqf(n,x,wsave)

******************************************************************

subroutine dcosqf computes the fast fourier transform of quarter
wave data. that is , dcosqf computes the coefficients in a cosine
series representation with only odd wave numbers. the transform
is defined below at output parameter x

dcosqf is the unnormalized inverse of dcosqb -- a call of dcosqf
followed by a call of dcosqb will multiply the input sequence x
by 4*n.

the array wsave which is used by subroutine dcosqf must be
initialized by calling subroutine dcosqi(n,wsave).


input parameters

n       the length of the array x to be transformed.  the method
        is most efficient when n is a product of small primes.

x       an array which contains the sequence to be transformed

wsave   a work array which must be dimensioned at least 3*n+15
        in the program that calls dcosqf. the wsave array must be
        initialized by calling subroutine dcosqi(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.

output parameters

x       for i=1,...,n

             x(i) = x(1) plus the sum from k=2 to k=n of

                2*x(k)*cos((2*i-1)*(k-1)*pi/(2*n))

             a call of dcosqf followed by a call of
             dcosqb will multiply the sequence x by 4*n.
             therefore dcosqb is the unnormalized inverse
             of dcosqf.

wsave   contains initialization calculations which must not
        be destroyed between calls of dcosqf or dcosqb.

******************************************************************

subroutine dcosqb(n,x,wsave)

******************************************************************

subroutine dcosqb computes the fast fourier transform of quarter
wave data. that is , dcosqb computes a sequence from its
representation in terms of a cosine series with odd wave numbers.
the transform is defined below at output parameter x.

dcosqb is the unnormalized inverse of dcosqf -- a call of dcosqb
followed by a call of dcosqf will multiply the input sequence x
by 4*n.

the array wsave which is used by subroutine dcosqb must be
initialized by calling subroutine dcosqi(n,wsave).


input parameters

n       the length of the array x to be transformed.  the method
        is most efficient when n is a product of small primes.

x       an array which contains the sequence to be transformed

wsave   a work array that must be dimensioned at least 3*n+15
        in the program that calls dcosqb. the wsave array must be
        initialized by calling subroutine dcosqi(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.

output parameters

x       for i=1,...,n

             x(i)= the sum from k=1 to k=n of

               4*x(k)*cos((2*k-1)*(i-1)*pi/(2*n))

             a call of dcosqb followed by a call of
             dcosqf will multiply the sequence x by 4*n.
             therefore dcosqf is the unnormalized inverse
             of dcosqb.

wsave   contains initialization calculations which must not
        be destroyed between calls of dcosqb or dcosqf.

******************************************************************

subroutine dcffti(n,wsave)

******************************************************************

subroutine dcffti initializes the array wsave which is used in
both dcfftf and dcfftb. the prime factorization of n together with
a tabulation of the trigonometric functions are computed and
stored in wsave.

input parameter

n       the length of the sequence to be transformed

output parameter

wsave   a work array which must be dimensioned at least 4*n+15
        the same work array can be used for both dcfftf and dcfftb
        as long as n remains unchanged. different wsave arrays
        are required for different values of n. the contents of
        wsave must not be changed between calls of dcfftf or dcfftb.

******************************************************************

subroutine dcfftf(n,c,wsave)

******************************************************************

subroutine dcfftf computes the forward complex discrete fourier
transform (the fourier analysis). equivalently , dcfftf computes
the fourier coefficients of a complex periodic sequence.
the transform is defined below at output parameter c.

the transform is not normalized. to obtain a normalized transform
the output must be divided by n. otherwise a call of dcfftf
followed by a call of dcfftb will multiply the sequence by n.

the array wsave which is used by subroutine dcfftf must be
initialized by calling subroutine dcffti(n,wsave).

input parameters


n      the length of the complex sequence c. the method is
       more efficient when n is the product of small primes. n

c      a complex array of length n which contains the sequence

wsave   a real work array which must be dimensioned at least 4n+15
        in the program that calls dcfftf. the wsave array must be
        initialized by calling subroutine dcffti(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.
        the same wsave array can be used by dcfftf and dcfftb.

output parameters

      for j=1,...,n

           c(j)=the sum from k=1,...,n of

                 c(k)*exp(-i*j*k*2*pi/n)

                       where i=sqrt(-1)

wsave   contains initialization calculations which must not be
        destroyed between calls of subroutine dcfftf or dcfftb

******************************************************************

subroutine dcfftb(n,c,wsave)

******************************************************************

subroutine dcfftb computes the backward complex discrete fourier
transform (the fourier synthesis). equivalently , dcfftb computes
a complex periodic sequence from its fourier coefficients.
the transform is defined below at output parameter c.

a call of dcfftf followed by a call of dcfftb will multiply the
sequence by n.

the array wsave which is used by subroutine dcfftb must be
initialized by calling subroutine dcffti(n,wsave).

input parameters


n      the length of the complex sequence c. the method is
       more efficient when n is the product of small primes.

c      a complex array of length n which contains the sequence

wsave   a real work array which must be dimensioned at least 4n+15
        in the program that calls dcfftb. the wsave array must be
        initialized by calling subroutine dcffti(n,wsave) and a
        different wsave array must be used for each different
        value of n. this initialization does not have to be
        repeated so long as n remains unchanged thus subsequent
        transforms can be obtained faster than the first.
        the same wsave array can be used by dcfftf and dcfftb.

output parameters

      for j=1,...,n

           c(j)=the sum from k=1,...,n of

                 c(k)*exp(i*j*k*2*pi/n)

                       where i=sqrt(-1)

wsave   contains initialization calculations which must not be
        destroyed between calls of subroutine dcfftf or dcfftb