DFFTPACK V1.0
*****************************************************************
A Double precision clone by Hugh C. Pumphrey of:
FFTPACK
version 4 april 1985
a package of fortran subprograms for the fast fourier
transform of periodic and other symmetric sequences
by
paul n swarztrauber
national center for atmospheric research boulder,colorado 80307
which is sponsored by the national science foundation
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
this package consists of programs which perform fast fourier
transforms for both double complex and (double precision) real
periodic sequences and certain other symmetric sequences that are
listed below.
1. dffti initialize dfftf and dfftb
2. dfftf forward transform of a real periodic sequence
3. dfftb backward transform of a real coefficient array
4. dzffti initialize dzfftf and dzfftb
5. dzfftf a simplified real periodic forward transform
6. dzfftb a simplified real periodic backward transform
7. dsinti initialize dsint
8. dsint sine transform of a real odd sequence
9. dcosti initialize dcost
10. dcost cosine transform of a real even sequence
11. dsinqi initialize dsinqf and dsinqb
12. dsinqf forward sine transform with odd wave numbers
13. dsinqb unnormalized inverse of dsinqf
14. dcosqi initialize dcosqf and dcosqb
15. dcosqf forward cosine transform with odd wave numbers
16. dcosqb unnormalized inverse of dcosqf
17. zffti initialize zfftf and zfftb
18. zfftf forward transform of a double complex periodic sequence
19. zfftb unnormalized inverse of zfftf
******************************************************************
subroutine dffti(n,wsave)
****************************************************************
subroutine dffti initializes the array wsave which is used in
both dfftf and dfftb. 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 dfftf and dfftb
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 dfftf or dfftb.
******************************************************************
subroutine dfftf(n,r,wsave)
******************************************************************
subroutine dfftf 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 dfftf. the wsave array must be
initialized by calling subroutine dffti(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 dfftf and dfftb.
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 dfftf
followed by a call of dfftb will multiply the input
sequence by n.
wsave contains results which must not be destroyed between
calls of dfftf or dfftb.
******************************************************************
subroutine dfftb(n,r,wsave)
******************************************************************
subroutine dfftb 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 dfftb. the wsave array must be
initialized by calling subroutine dffti(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 dfftf and dfftb.
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 dfftf
followed by a call of dfftb will multiply the input
sequence by n.
wsave contains results which must not be destroyed between
calls of dfftb or dfftf.
******************************************************************
subroutine dzffti(n,wsave)
******************************************************************
subroutine dzffti initializes the array wsave which is used in
both dzfftf and dzfftb. 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 3*n+15.
the same work array can be used for both dzfftf and dzfftb
as long as n remains unchanged. different wsave arrays
are required for different values of n.
******************************************************************
subroutine dzfftf(n,r,azero,a,b,wsave)
******************************************************************
subroutine dzfftf computes the fourier coefficients of a real
perodic sequence (fourier analysis). the transform is defined
below at output parameters azero,a and b. dzfftf is a simplified
but slower version of dfftf.
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 which must be dimensioned at least 3*n+15.
in the program that calls dzfftf. the wsave array must be
initialized by calling subroutine dzffti(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 dzfftf and dzfftb.
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 dzfftb(n,r,azero,a,b,wsave)
******************************************************************
subroutine dzfftb computes a real perodic sequence from its
fourier coefficients (fourier synthesis). the transform is
defined below at output parameter r. dzfftb is a simplified
but slower version of dfftb.
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 which must be dimensioned at least 3*n+15.
in the program that calls dzfftb. the wsave array must be
initialized by calling subroutine dzffti(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 dzfftf and dzfftb.
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
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 since 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 since 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 since 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
cosqb 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 since 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 zffti(n,wsave)
******************************************************************
subroutine zffti initializes the array wsave which is used in
both zfftf and zfftb. 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 zfftf and zfftb
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 zfftf or zfftb.
******************************************************************
subroutine zfftf(n,c,wsave)
******************************************************************
subroutine zfftf computes the forward complex discrete fourier
transform (the fourier analysis). equivalently , zfftf 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 zfftf
followed by a call of zfftb will multiply the sequence by n.
the array wsave which is used by subroutine zfftf must be
initialized by calling subroutine zffti(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 zfftf. the wsave array must be
initialized by calling subroutine zffti(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 zfftf and zfftb.
output parameters
c for j=1,...,n
c(j)=the sum from k=1,...,n of
c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
where i=sqrt(-1)
wsave contains initialization calculations which must not be
destroyed between calls of subroutine zfftf or zfftb
******************************************************************
subroutine zfftb(n,c,wsave)
******************************************************************
subroutine zfftb computes the backward complex discrete fourier
transform (the fourier synthesis). equivalently , zfftb computes
a complex periodic sequence from its fourier coefficients.
the transform is defined below at output parameter c.
a call of zfftf followed by a call of zfftb will multiply the
sequence by n.
the array wsave which is used by subroutine zfftb must be
initialized by calling subroutine zffti(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 zfftb. the wsave array must be
initialized by calling subroutine zffti(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 zfftf and zfftb.
output parameters
c for j=1,...,n
c(j)=the sum from k=1,...,n of
c(k)*exp(i*(j-1)*(k-1)*2*pi/n)
where i=sqrt(-1)
wsave contains initialization calculations which must not be
destroyed between calls of subroutine zfftf or zfftb
["send index for vfftpk" describes a vectorized version of fftpack]