Currently, the SCSL is available only on SGI machines.
In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. Suppose that the arrays are dimensioned as follows:
COMPLEX X(0:N-1), Y(0:N-1)
The output array is the FFT of the input array, using the following formula for the FFT:
n-1
Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ] for k = 0, ..., n-1
j=0
where:
w = exp(2*pi*i/n),
i = + sqrt(-1),
pi = 3.14159...,
isign = +1 or -1
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make this routine compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any particular values of isign and scale, the mathematical inverse function is computed by taking the FFT with -isign and 1/(n*scale). In particular, if you use isign = +1 and scale = 1.0 you can compute the inverse FFT by using isign = -1 and scale = 1.0/n.
The output array may be the same as the input array.
COMPLEX X(0:N-1)
COMPLEX Y(0:N-1)
However, if you prefer to use the more customary FORTRAN style with subscripts starting at 1 you do not have to change the calling sequence, as in the following (assuming N > 0):
COMPLEX X(N)
COMPLEX Y(N)
Example 1: Initialize the complex array table in preparation for doing an FFT of size 1024. Only the isign, n, and table arguments are used in this case. You can use dummy arguments or zeros for the other arguments in the subroutine call.
REAL TABLE(30 + 2048)
CALL CCFFT(0, 1024, 0.0, DUMMY, DUMMY, TABLE, DUMMY, 0)
Example 2: x and y are complex arrays of dimension (0:1023). Take the FFT of x and store the results in y. Before taking the FFT, initialize the table array, as in example 1.
COMPLEX X(0:1023), Y(0:1023)
REAL TABLE(30 + 2048)
REAL WORK(2048)
CALL CCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
CALL CCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
Example 3: Using the same x and y as in example 2, take the inverse FFT of y and store it back in x. The scale factor 1/1024 is used. Assume that the table array is already initialized.
CALL CCFFT(-1, 1024, 1.0/1024.0, Y, X, TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change was needed in the subroutine calls.
COMPLEX X(1024), Y(1024)
CALL CCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
CALL CCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
Example 5: Do the same computation as in example 4, but put the output back in array x to save storage space. Assume that table is already initialized.
COMPLEX X(1024)
CALL CCFFT(1, 1024, 1.0, X, X, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
Input array of values to be transformed.
Table of factors and trigonometric functions.
If isign = 0, the routine initializes table (table is output only).
If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0 (table is input only).
Work array. This is a scratch array used for intermediate calculations. Its address space must be different address space from that of the input and output arrays.
scfft/dzfft computes the FFT of the real array x, and it stores the results in the complex array y. csfft/zdfft computes the corresponding inverse complex-to-real transform.
It is customary in FFT applications to use zero-based subscripts; the formulas are simpler that way. For these routines, suppose that the arrays are dimensioned as follows:
REAL X(0:n-1)
COMPLEX Y(0:n/2)
Then the output array is the FFT of the input array, using the following formula for the FFT:
n-1
Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ] for k = 0, ..., n/2
j=0
where:
w = exp(2*pi*i/n),
i = + sqrt(-1),
pi = 3.14159...,
isign = +1 or -1.
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make these routines compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you call scfft with any particular values of isign and scale, the mathematical inverse function is computed by calling csfft with -isign and 1/(n*scale). In particular, if you use isign = +1 and scale = 1.0 in scfft for the forward FFT, you can compute the inverse FFT by using ccfft with isign = -1 and scale = 1.0/n.
The mathematical definition of the Fourier transform takes a sequence of n complex values and transforms it to another sequence of n complex values. A complex-to-complex FFT routine, such as ccfft, will take n complex input values, and produce n complex output values. In fact, one easy way to compute a real-to-complex FFT is to store the input data in a complex array, then call routine ccfft to compute the FFT. You get the same answer when using the scfft routine.
The reason for having a separate real-to-complex FFT routine is efficiency. Because the input data is real, you can make use of this fact to save almost half of the computational work. The theory of Fourier transforms tells us that for real input data, you have to compute only the first n/2 + 1 complex output values, because the remaining values can be computed from the first half of the values by the simple formula:
Y(k) = conjg(Y(n-k)) for n/2 <= k <= n-1
where the notation conjgY represents the complex conjugate of y.
In fact, in many applications, the second half of the complex output data is never explicitly computed or stored. Likewise, as explained later, only the first half of the complex data has to be supplied for the complex-to-real FFT.
Another implication of FFT theory is that, for real input data, the first output value, Y(0), will always be a real number; therefore, the imaginary part will always be 0. If n is an even number, Y(n/2) will also be real and thus, have zero imaginary parts.
Generally, the FFT transforms a complex sequence into a complex sequence. However, in a certain application we may know the output sequence is real. Often, this is the case because the complex input sequence was the transform of a real sequence. In this case, you can save about half of the computational work.
According to the theory of Fourier transforms, for the output sequence, Y, to be a real sequence, the following identity on the input sequence, X, must be true:
X(k) = conjg(X(n-k)) for n/2 <= k <= n-1
And, in fact, the input values X(k) for k > n/2 need not be supplied; they can be inferred from the first half of the input.
Thus, in the complex-to-real routine, csfft, the arrays can be dimensioned as follows:
COMPLEX X(0:n/2)
REAL Y(0:n-1)
There are n/2 + 1 complex input values and n real output values. Even though only n/2 + 1 input values are supplied, the size of the transform is still n in this case, because implicitly you are using the FFT formula for a sequence of length n.
Another implication of the theory is that X(0) must be a real number (that is, it must have zero imaginary part). Also, if n is even, X(n/2) must also be real. Routine CSFFT assumes that these values are real; if you specify a nonzero imaginary part, it is ignored.
REAL X(0:n-1)
COMPLEX Y(0:n/2)
No change is needed in the calling sequence; however, if you prefer you can use the more customary Fortran style with subscripts starting at 1, as in the following:
REAL X(n)
COMPLEX Y(n/2 + 1)
Example 1: Initialize the complex array TABLE in preparation for doing an FFT of size 1024. In this case only the arguments isign, n, and table are used. You can use dummy arguments or zeros for the other arguments in the subroutine call.
REAL TABLE(15 + 1024)
CALL SCFFT(0, 1024, 0.0, DUMMY, DUMMY, TABLE, DUMMY, 0)
Example 2: X is a real array of dimension (0:1023), and Y is a complex array of dimension (0:512). Take the FFT of X and store the results in Y. Before taking the FFT, initialize the TABLE array, as in example 1.
REAL X(0:1023)
COMPLEX Y(0:512)
REAL TABLE(15 + 1024)
REAL WORK(1024)
CALL SCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1/1024 is used. Assume that the TABLE array is initialized already.
CALL CSFFT(-1, 1024, 1.0/1024.0, Y, X, TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls are not changed.
REAL X(1024)
COMPLEX Y(513)
CALL SCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. Assume that the TABLE array is initialized already.
REAL X(1024)
COMPLEX Y(513)
EQUIVALENCE ( X(1), Y(1) )
CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0. If isys(0) > 0, isys(0) gives the upper bound of the isys array; that is, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).
The output array, y, is the FFT of the the input array, x, computed according to the preceding formula. The output array may be equivalenced to the input array in the calling program. Be careful when dimensioning the arrays, in this case, to allow for the fact that the complex array contains two (real) words more than the real array.
Table of factors and trigonometric functions.
If isign = 0, the table array is initialized to contain trigonometric tables needed to compute an FFT of size n.
If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0.
Work array used for intermediate calculations. Its address space must be different from that of the input and output arrays.
ccfftm/zzfftm computes the FFT of each column of the complex matrix x, and stores the results in the columns of complex matrix y.
Suppose the arrays are dimensioned as follows:
COMPLEX X(0:ldx-1, 0:lot-1)
COMPLEX Y(0:ldy-1, 0:lot-1)
where ldx >= n, ldy >= n.
Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:
n-1
Y(k, L) = scale * Sum [ X(j)*w**(isign*j*k) ]
j=0
for k = 0, ..., n-1
L = 0, ..., lot-1
where:
w = exp(2*pi*i/n),
i = + sqrt(-1),
pi = 3.14159...,
isign = +1 or -1
lot = the number of columns to transform
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make this routine compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any particular values of isign and scale, the mathematical inverse function is computed by taking the FFT with -isign and 1/(n * scale). In particular, if you use isign = +1 and scale = 1.0 for the forward FFT, you can compute the inverse FFT by using the following: isign = -1 and scale = 1.0/n.
This section contains information about the algorithm for these routines, the initialization of the table array, the declaration of dimensions for x and y arrays, some performance tips, and some implementation-dependent details.
COMPLEX X(0:ldx-1, 0:lot-1)
COMPLEX Y(0:ldy-1, 0:lot-1)
The calling sequence does not have to change, however, if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:
COMPLEX X(ldx, lot)
COMPLEX Y(ldy, lot)
REAL TABLE(30 + 256)
CALL CCFFTM(0, 128, 0, 0., DUMMY, 1, DUMMY, 1, TABLE, DUMMY, 0)
Example 2: X and Y are complex arrays of dimension (0:128) by (0:55). The first 128 elements of each column contain data. For performance reasons, the extra element forces the leading dimension to be an odd number. Take the FFT of the first 50 columns of X and store the results in the first 50 columns of Y. Before taking the FFT, initialize the TABLE array, as in example 1.
COMPLEX X(0:128, 0:55)
COMPLEX Y(0:128, 0:55)
REAL TABLE(30 + 256)
REAL WORK(256)
...
CALL CCFFTM(0, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
CALL CCFFTM(1, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1/128 is used. Assume that the TABLE array is already initialized.
CALL CCFFTM(-1, 128, 50, 1./128., Y, 129, X, 129, TABLE,WORK,0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls are not changed.
COMPLEX X(129, 55)
COMPLEX Y(129, 55)
...
CALL CCFFTM(0, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
CALL CCFFTM(1, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but put the output back in array X to save storage space. Assume that the TABLE array is already initialized.
COMPLEX X(129, 55)
...
CALL CCFFTM(1, 128, 50, 1.0, X, 129, X, 129, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0.
If isys(0) > 0, isys(0) gives the upper bound of the isys array. Therefore, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).
The output array may be the same as the input array. In that case, the transform is done in place. The input array is overwritten with the transformed values. In this case, it is necessary that ldx = ldy.
If isign = 0, the routine initializes table (table is output only).
If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0 (table is input only).
scfftm/dzfftm computes the FFT of each column of the real matrix X, and it stores the results in the corresponding column of the complex matrix Y. csfftm/zdfftm computes the corresponding inverse transforms.
In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. First, the function of scfftm is described. Suppose that the arrays are dimensioned as follows:
REAL X(0:ldx-1, 0:lot-1)
COMPLEX Y(0:ldy-1, 0:lot-1)
where ldx >= n, ldy >= n/2 + 1.
Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:
n-1
Y(k, L) = scale * Sum [ X(j, L)*w**(isign*j*k) ]
j=0
for k = 0, ..., n/2
L = 0, ..., lot-1 where:
w = exp(2*pi*i/n),
i = + sqrt(-1)
pi = 3.14159...,
isign = +1 or -1,
lot = the number of columns to transform
Different authors use different conventions for which transform (isign = +1 or isign = -1) is used in the real-to-complex case, and what the scale factor should be. Some adopt the convention that isign = 1 for the real-to-complex transform, and isign = -1 for the complex-to-real inverse. Others use the opposite convention. You can make these routines compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you use scfftm to take the real-to-complex FFT, using any particular values of isign and scale, the mathematical inverse function is computed by using csfftm with -isign and 1/ (n*scale). In particular, if you call scfftm with isign = +1 and scale = 1.0, you can use csfftm to compute the inverse complex-to-real FFT by using isign = -1 and scale = 1.0/n.
The mathematical definition of the Fourier transform takes a sequence of n complex values and transforms it to another sequence of n complex values. A complex-to-complex FFT routine, such as ccfftm, will take n complex input values and produce n complex output values. In fact, one easy way to compute a real-to-complex FFT is to store the input data x in a complex array, then call routine ccfftm to compute the FFT. You get the same answer when using the scfftm routine.
A separate real-to-complex FFT routine is more efficient than the equivalent complex-to-complex routine. Because the input data is real, you can make use of this fact to save almost half of the computational work. According to the theory of Fourier transforms, for real input data, you have to compute only the first n/2 + 1 complex output values in each column, because the second half of the FFT values in each column can be computed from the first half of the values by the simple formula:
Y = conjgY for n/2 <= k <= n-1
k,L n-k, L
where the notation conjg(z) represents the complex conjugate of z.
In fact, in many applications, the second half of the complex output data is never explicitly computed or stored. Likewise, you must supply only the first half of the complex data in each column has to be supplied for the complex-to-real FFT.
Another implication of FFT theory is that for real input data, the first output value in each column, Y(0, L), will always be a real number; therefore, the imaginary part will always be 0. If n is an even number, Y(n/2, L) will also be real and have 0 imaginary parts.
In general, the FFT transforms a complex sequence into a complex sequence; however, in a certain application you may know the output sequence is real, perhaps because the complex input sequence was the transform of a real sequence. In this case, you can save about half of the computational work.
According to the theory of Fourier transforms, for the output sequence, Y, to be a real sequence, the following identity on the input sequence, X, must be true:
X = conjgX for n/2 <= k <= n-1
k,L n-k,L
And, in fact, the following input values
X for k > n/2
k,L
do not have to be supplied, because they can be inferred from the
first half of the input.
Thus, in the complex-to-real routine, CSFFTM, the arrays can be dimensioned as follows:
COMPLEX X(0:ldx-1, 0:lot-1)
REAL Y(0:ldy-1, 0:lot-1)
where ldx >= n/2 + 1, ldy >= n.
In each column, there are (n/2) + 1 complex input values and n real output values. Even though only (n/2) + 1 input values are supplied, the size of the transform is still n in this case, because implicitly the FFT formula for a sequence of length n is used.
Another implication of the theory is that X(0, L) must be a real number (that is, must have zero imaginary part). If n is an even number, X(n/2, L) must also be real. Routine CSFFTM assumes that each of these values is real; if a nonzero imaginary part is given, it is ignored.
REAL X(0:ldx-1, 0:lot-1)
COMPLEX Y(0:ldy-1, 0:lot-1)
No change is made in the calling sequence, however, if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:
REAL X(ldx, lot)
COMPLEX Y(ldy, lot)
Example 1: Initialize the complex array TABLE in preparation for doing an FFT of size 128. In this case only the isign, n, and table arguments are used; you may use dummy arguments or zeros for the other arguments in the subroutine call.
REAL TABLE(15 + 128)
CALL SCFFTM(0, 128, 1, 0.0, DUMMY, 1, DUMMY, 1,
& TABLE, DUMMY, 0)
Example 2: X is a real array of dimension (0:128, 0:55), and Y is a complex array of dimension (0:64, 0:55). The first 128 elements in each column of X contain data; the extra element forces an odd leading dimension. Take the FFT of the first 50 columns of X and store the results in the first 50 columns of Y. Before taking the FFT, initialize the TABLE array, as in example 1.
REAL X(0:128, 0:55)
COMPLEX Y(0:64, 0:55)
REAL TABLE(15 + 128)
REAL WORK((128)
...
CALL SCFFTM(0, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
CALL SCFFTM(1, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1/128 is used. Assume that the TABLE array is initialized already.
CALL CSFFTM(-1, 128, 50, 1.0/128.0, Y, 65, X, 129,
& TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change is made in the subroutine calls.
REAL X(129, 56)
COMPLEX Y(65, 56)
...
CALL SCFFTM(0, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
CALL SCFFTM(1, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. In this case, a row must be added to X, because it is equivalenced to a complex array. The leading dimension of X is two times an odd number; therefore, memory bank conflicts are minimal. Assume that TABLE is initialized already.
REAL X(130, 56)
COMPLEX Y(65, 56)
EQUIVALENCE ( X(1, 1), Y(1, 1) )
...
CALL SCFFTM(1, 128, 50, 1.0, X, 130, Y, 65, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0.
If isys(0) > 0, isys(0) gives the upper bound of the isys array. Therefore, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).
Each column of the output array, y, is the FFT of the corresponding column of the input array, x, computed according to the preceding formula. The output array may be equivalenced to the input array. In that case, the transform is done in place and the input array is overwritten with the transformed values. In this case, the following conditions on the leading dimensions must hold:
scfftm, dzfftm: ldx = 2ldy. csfftm, zdfftm: ldy = 2ldx.
If isign = 0, table is initialized to contain trigonometric tables needed to compute an FFT of length n.
These routines compute the two-dimensional complex Fast Fourier Transform (FFT) of the complex matrix x, and store the results in the complex matrix y. ccfft2d does the complex-to-complex transform and zzfft does the same for double precision arrays.
In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. Suppose that the arrays are dimensioned as follows:
COMPLEX X(0:n1-1, 0:n2-1)
COMPLEX Y(0:n1-1, 0:n2-1)
These routines compute the formula:
n2-1 n1-1
Y(k1, k2) = scale * Sum Sum [ X(j1, j2)*w1**(j1*k1)*w2**(j2*k2) ]
j2=0 j1=0
for k1 = 0, ..., n1-1
k2 = 0, ..., n2-1
where:
w1 = exp(isign*2*pi*i/n1)
w2 = exp(isign*2*pi*i/n2)
i = + sqrt(-1)
pi = 3.14159...,
isign = +1 or -1
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make this routine compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any particular values of isign and scale, the mathematical inverse function is computed by taking the FFT with -isign and 1/(n1*n2*scale). In particular, if you use isign = +1 and scale = 1.0 for the forward FFT, you can compute the inverse FFT by using isign = -1 and scale = 1.0/(n1*n2).
COMPLEX X(0:ldx-1, 0:n2-1)
COMPLEX Y(0:ldy-1, 0:n2-1)
However, the calling sequence does not change if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:
COMPLEX X(ldx, n2)
COMPLEX Y(ldy, n2)
Example 1: Initialize the TABLE array in preparation for doing a two-dimensional FFT of size 128 by 256. In this case only the isign, n1, n2, and table arguments are used; you can use dummy arguments or zeros for other arguments.
REAL TABLE ((30 + 256) + (30 + 512))
CALL CCFFT2D (0, 128, 256, 0.0, DUMMY, 1, DUMMY, 1,
& TABLE, DUMMY, 0)
Example 2: X and Y are complex arrays of dimension (0:128, 0:255). The first 128 elements of each column contain data. For performance reasons, the extra element forces the leading dimension to be an odd number. Take the two-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
COMPLEX X(0:128, 0:255)
COMPLEX Y(0:128, 0:255)
REAL TABLE((30 + 256) + (30 + 512))
REAL WORK(2*128*256)
...
CALL CCFFT2D(0, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
CALL CCFFT2D(1, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1/(128*256) is used. Assume that the TABLE array is already initialized.
CALL CCFFT2D(-1, 128, 256, 1.0/(128.0*256.0), Y, 129,
& X, 129, TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls are not changed.
COMPLEX X(129, 256)
COMPLEX Y(129, 256)
...
CALL CCFFT2D(0, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
CALL CCFFT2D(1, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but put the output back in array X to save storage space. Assume that the TABLE array is already initialized.
COMPLEX X(129, 256)
...
CALL CCFFT2D(1, 128, 256, 1.0, X, 129, X, 129, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n1, n2, table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
The number of rows in the y array, as it was declared in the calling program (the leading dimension of y). ldy >= MAX(n1, 1).
Table of factors and trigonometric functions.
If isign = 0, the routine initializes table (table is output only).
If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0 (table is input only).
Work array. This is a scratch array used for intermediate calculations. Its address space must be different from that of the input and output arrays.
scfft2d/dzfft2d computes the two-dimensional Fast Fourier Transform (FFT) of the real matrix x, and it stores the results in the complex matrix y. csfft2d/zdfft2d computes the corresponding inverse transform.
In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. First the function of scfft2d is described. Suppose the arrays are dimensioned as follows:
REAL X(0:ldx-1, 0:n2-1)
COMPLEX Y(0:ldy-1, 0:n2-1)
where ldx >= n1 ldy >= (n1/2) + 1.
scfft2d computes the formula:
n2-1 n1-1
Y(k1, k2) = scale * Sum Sum [ X(j1, j2)*w1**(j1*k1)*w2**(j2*k2) ]
j2=0 j1=0
for k1 = 0, ..., n1/2 + 1
k2 = 0, ..., n2-1
where:
w1 = exp(isign*2*pi*i/n1)
w2 = exp(isign*2*pi*i/n2)
i = + sqrt(-1)
pi = 3.14159...,
isign = +1 or -1
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make these routines compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any particular values of isign and scale, the mathematical inverse function is computed by taking the FFT with -isign and 1/(n1 * n2 * scale). In particular, if you use isign = +1 and scale = 1.0 for the forward FFT, you can compute the inverse FFT by using isign = -1 and scale = 1.0/(n1 . n2).
scfft2d is very similar in function to ccfft2d, but it takes the real-to-complex transform in the first dimension, followed by the complex-to-complex transform in the second dimension.
csfft2d does the reverse. It takes the complex-to-complex FFT in the second dimension, followed by the complex-to-real FFT in the first dimension.
See the scfft man page for more information about real-to-complex and complex-to-real FFTs. The two-dimensional analog of the conjugate formula is as follows:
Y = conjg Y
k , k n1 - k , n2 - k
1 2 1 2
for n1/2 < k <= n1 - 1
1
0 <= k <= n2 - 1
2
where the notation conjg(z) represents the complex conjugate of z.
Thus, you have to compute only (slightly more than) half of the output values, namely:
Y for 0 <= k <= n1/2 0 <= k <= n2 - 1
k , k 1 2
1 2
csfft2d uses a routine similar to ccfftm to do a complex-to-complex FFT on the rows, then uses a routine similar to csfftm to do a complex-to-real FFT on the columns.
REAL X(0:ldx-1, 0:n2-1)
COMPLEX Y(0:ldy-1, 0:n2-1)
No change is made in the calling sequence, however, if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:
REAL X(ldx, n2)
COMPLEX Y(ldy, n2)
Example 1: Initialize the TABLE array in preparation for doing a two-dimensional FFT of size 128 by 256. In this case, only the isign, n1, n2, and table arguments are used; you can use dummy arguments or zeros for other arguments.
REAL TABLE ((15 + 128) + 2(15 + 256))
CALL SCFFT2D (0, 128, 256, 0.0, DUMMY, 1, DUMMY, 1,
& TABLE, DUMMY, 0)
Example 2: X is a real array of size (0:128, 0: 255), and Y is a complex array of dimension (0:64, 0:255). The first 128 elements of each column of X contain data; for performance reasons, the extra element forces the leading dimension to be an odd number. Take the two-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
REAL X(0:128, 0:255)
COMPLEX Y(0:64, 0:255)
REAL TABLE ((15 + 128) + 2(15 + 256))
REAL WORK(128*256)
...
CALL SCFFT2D(0, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
CALL SCFFT2D(1, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1/(128*256) is used. Assume that the TABLE array is initialized already.
CALL CSFFT2D(-1, 128, 256, 1.0/(128.0*256.0), Y, 65,
& X, 130, TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change is needed in the subroutine calls.
REAL X(129, 256)
COMPLEX Y(65, 256)
...
CALL SCFFT2D(0, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
CALL SCFFT2D(1, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. In this case, a row must be added to X, because it is equivalenced to a complex array. Assume that TABLE is already initialized.
REAL X(130, 256)
COMPLEX Y(65, 256)
EQUIVALENCE ( X(1, 1), Y(1, 1) )
...
CALL SCFFT2D(1, 128, 256, 1.0, X, 130, Y, 65, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
Array of values to be transformed.
The number of rows in the y array, as it was declared in the calling program (the leading dimension of y).
scfft2d, dzfft2d: ldy >= MAX(n1/2 + 1, 1). csfft2d, zdfft2d: ldy >= MAX(n1 + 2, 1).
In the complex-to-real routine, two extra elements are in the first dimension (ldy >= n1 + 2, rather than just ldy >= n1). These elements are needed for intermediate storage during the computation. On exit, their value is undefined.
If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0.
If isys(0) > 0, isys(0) gives the upper bound of the isys array. Therefore, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).
Output array of transformed values. The output array can be the same as the input array, in which case, the transform is done in place and the input array is overwritten with the transformed values. In this case, it is necessary that the following equalities hold:
scfft2d, dzfft2d: ldx = 2 * ldy. csfft2d, zdfft2d: ldy = 2 * ldx.
Table of factors and trigonometric functions.
If isign = 0, the routine initializes table (table is output only).
If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0 (table is input only).
Work array. This is a scratch array used for intermediate calculations. Its address space must be different from that of the input and output arrays.
These routines compute the three-dimensional complex FFT of the complex matrix X, and store the results in the complex matrix Y.
In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. So suppose the arrays are dimensioned as follows:
COMPLEX X(0:n1-1, 0:n2-1, 0:n3-1)
COMPLEX Y(0:n1-1, 0:n2-1, 0:n3-1)
These routines compute the formula:
Y(k1,k2,k3) =
n1-1 n2-1 n3-1
scale * Sum Sum Sum [X(j1,j2,j3)*w1**(j1*k1)*w2**(j2*k2)*w3**(j3*k3)]
j1=0 j2=0 j3=0
for k1 = 0, ..., n1 - 1,
k2 = 0, ..., n2 - 1,
k3 = 0, ..., n3 - 1,
where:
w1 = exp(isign*2*pi*i/n1),
w2 = exp(isign*2*pi*i/n2),
w3 = exp(isign*2*pi*i/n3),
i = + sqrt(-1)
pi = 3.14159...
isign = +1 or -1
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make this routine compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any particular values of isign and scale, the mathematical inverse function is computed by taking the FFT with -isign and 1/(n1 * n2 * n3 * scale). In particular, if you use isign = +1 and scale = 1.0 for the forward FFT, you can compute the inverse FFT by using isign = -1 and scale = 1/(n1 . n2 . n3).
Example 1: Initialize the TABLE array in preparation for doing a three-dimensional FFT of size 128 by 128 by 128. In this case, only the isign, n1, n2, n3, and table arguments are used; you can use dummy arguments or zeros for other arguments.
REAL TABLE ((30 + 256) + (30 + 256) + (30 + 256))
CALL CCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
& TABLE, DUMMY, 0)
Example 2: X and Y are complex arrays of dimension (0:128, 0:128, 0:128). The first 128 elements of each dimension contain data; for performance reasons, the extra element forces the leading dimensions to be odd numbers. Take the three-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
COMPLEX X(0:128, 0:128, 0:128)
COMPLEX Y(0:128, 0:128, 0:128)
REAL TABLE ((30+256) + (30 + 256) + (30 + 256))
REAL WORK 2(128*128*128)
...
CALL CCFFT3D(0, 128, 128, 128, 1.0, DUMMY, 1, 1,
& DUMMY, 1, 1, TABLE, WORK, 0)
CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
& Y, 129, 129, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1.0/(128.0**3) is used. Assume that the TABLE array is already initialized.
CALL CCFFT3D(-1, 128, 128, 128, 1.0/(128.0**3), Y, 129, 129,
& X, 129, 129, TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls do not change.
COMPLEX X(129, 129, 129)
COMPLEX Y(129, 129, 129)
...
CALL CCFFT3D(0, 128, 128, 128, 1.0, DUMMY, 1, 1,
& DUMMY, 1, 1, TABLE, WORK, 0)
CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
& Y, 129, 129, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but put the output back in the array X to save storage space. Assume that the TABLE array is already initialized.
COMPLEX X(129, 129, 129)
...
CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
& X, 129, 129, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n1, n2, n3, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
Each element of the output array is multiplied by scale after taking the Fourier transform, as defined previously.
Input array of values to be transformed.
isys = 0 or 1 depending on the amount of workspace the user can provide to the routine.
Output array of transformed values. The output array may be the same as the input array, in which case, the transform is done in place; that is, the input array is overwritten with the transformed values. In this case, it is necessary that ldx = ldy, and ldx2 = ldy2.
Table of factors and trigonometric functions. If isign = 0, the routine initializes table (table is output only). If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0 (table is input only).
Work array. This is a scratch array used for intermediate calculations. Its address space must be different from that of the input and output arrays.
These are C++ wrapper functions for the 3D real-to-complex and complex-to-real transform routines in the SGI/Cray Scientific Library (SCSL). The purpose of these definitions is to overload the functions so that C++ users can access the functions in SCSL with identical function names.
scfft3d/dzfft3d computes the three-dimensional Fast Fourier Transform (FFT) of the real matrix X, and it stores the results in the complex matrix Y. csfft3d/zdfft3d computes the corresponding inverse transform.
In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. First, the function of SCFFT3D is described. Suppose the arrays are dimensioned as follows:
REAL X(0:ldx-1, 0:ldx2-1, 0:n3-1)
COMPLEX Y(0:ldy-1, 0:ldy2-1, 0:n3-1)
scfft3d computes the formula:
Y(k1,k2,k3) =
n1-1 n2-1 n3-1
scale * Sum Sum Sum [X(j1,j2,j3)*w1**(j1*k1)*w2**(j2*k2)*w3**(j3*k3)]
j1=0 j2=0 j3=0
for k1 = 0, ..., n1/2,
k2 = 0, ..., n2 - 1,
k3 = 0, ..., n3 - 1,
where:
w1 = exp(isign*2*pi*i/n1),
w2 = exp(isign*2*pi*i/n2),
w3 = exp(isign*2*pi*i/n3),
i = + sqrt(-1)
pi = 3.14159...
isign = +1 or -1
Different authors use different conventions for which of the transforms, isign = +1 or isign = -1, is the forward or inverse transform, and what the scale factor should be in either case. You can make these routines compute any of the various possible definitions, however, by choosing the appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any particular values of isign and scale, the mathematical inverse function is computed by taking the FFT with -isign and 1/(n1 * n2 * n3 * scale). In particular, if you use isign = +1 and scale = 1.0 for the forward FFT, you can compute the inverse FFT by isign = -1 and
scale = 1.0/(n1*n2*n3).
scfft3d is very similar in function to ccfft3d, but it takes the real-to-complex transform in the first dimension, followed by the complex-to-complex transform in the second and third dimensions.
csfft3d does the reverse. It takes the complex-to-complex FFT in the third and second dimensions, followed by the complex-to-real FFT in the first dimension.
See the scfftm man page for more information about real-to-complex and complex-to-real FFTs. The three dimensional analog of the conjugate formula is as follows:
Y = conjg Y
k ,k ,k n1 - k , n2 - k , n3 - k
1 2 3 1 2 3
for n1/2 < k <= n1 - 1
1
0 <= k <= n2 - 1
2
0 <= k <= n3 - 1
3
where the notation conjg(z) represents the complex conjugate of z.
Thus, you have to compute only (slightly more than) half out the output values, namely:
Y
k ,k ,k
1 2 3
for 0 <= k <= n1/2
1
0 <= k <= n2 - 1
2
0 <= k <= n3 - 1
The following examples are for Origin series only.
Example 1: Initialize the TABLE array in preparation for doing a three-dimensional FFT of size 128 by 128 by 128. In this case only the isign, n1, n2, n3, and table arguments are used; you can use dummy arguments or zeros for other arguments.
REAL TABLE ((15 + 128) + 2(15+128) + 2( 15 + 128))
CALL SCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
& TABLE, DUMMY, 0)
Example 2: X is a real array of size (0:128, 0:128, 0:128). The first 128 elements of each dimension contain data; for performance reasons, the extra element forces the leading dimensions to be odd numbers. Y is a complex array of dimension (0:64, 0:128, 0:128). Take the three-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
REAL X(0:128, 0:128, 0:128)
COMPLEX Y(0:64, 0:128, 0:128)
REAL TABLE ((15+128) + 2(15 + 128) + 2(15 + 128))
REAL WORK(128*128*128)
...
CALL SCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
& Y, 65, 129, TABLE, WORK, 0)
CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
& Y, 65, 129, TABLE, WORK, 0)
Example 3: With X and Y as in example 2, take the inverse FFT of Y and store it back in X. The scale factor 1/(128**3) is used. Assume that the TABLE array is initialized already.
CALL CSFFT3D(-1, 128, 128, 128, 1.0/128.0**3, Y, 65, 129,
& X, 130, 129, TABLE, WORK, 0)
Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change is made in the subroutine calls.
REAL X(129, 129, 129)
COMPLEX Y(65, 129, 129)
REAL TABLE ((15+128) + 2(15 + 128) + 2(15 + 128))
REAL WORK(128*128*128)
...
CALL SCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
& Y, 65, 129, TABLE, WORK, 0)
CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
& X, 129, 129, TABLE, WORK, 0)
Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. Assume that the TABLE array is initialized already.
REAL X(130, 129, 129)
COMPLEX Y(65, 129, 129)
EQUIVALENCE (X(1, 1, 1), Y(1, 1, 1))
...
CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 130, 129,
& Y, 65, 129, TABLE, WORK, 0)
Input parameters:
If isign = 0, the routine initializes the table array and returns. In this case, the only arguments used or checked are isign, n1, n2, n3, and table.
If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.
Array of values to be transformed.
scfft3d, dzfft3d: ldx >= MAX(n1, 1). csfft3d, zdfft3d: ldx >= MAX(n1/2 + 1, 1).
scfft3d, dzfft3d: ldy >= MAX(n1/2 + 1, 1). csfft3d, zdfft3d: ldy >= MAX(n1 + 2, 1).
In the complex-to-real routine, two extra elements are in the first dimension (that is, ldy >= n1 + 2, rather than just ldy >= n1). These elements are needed for intermediate storage during the computation. On exit, their value is undefined.
isys = 0 or 1 depending on the amount of workspace the user can provide to the routine.
Output array of transformed values. The output array can be the same as the input array, in which case, the transform is done in place; that is, the input array is overwritten with the transformed values. In this case, it is necessary that the following equalities hold:
scfft3d, dzfft3d: ldx = 2 * ldy, and ldx2 = ldy2. csfft3d, zdfft3d: ldy = 2 * ldx, and ldx2 = ldy2.
Table of factors and trigonometric functions.
This array must be initialized by a call to scfft3d or csfft3d with isign = 0.
If isign = 0, table is initialized to contain trigonometric tables needed to compute a three-dimensional FFT of size n1 by n2 by n3. If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0.
Work array. This is a scratch array used for intermediate calculations. Its address space must be different from that of the input and output arrays.