|
casacore
|
#include <SCSL.h>
Static Public Member Functions | |
| static void | ccfft (Int isign, Int n, Float scale, Complex *x, Complex *y, Float *table, Float *work, Int isys) |
| These routines compute the Fast Fourier Transform (FFT) of the complex vector x, and store the result in vector y. More... | |
| static void | ccfft (Int isign, Int n, Double scale, DComplex *x, DComplex *y, Double *table, Double *work, Int isys) |
| static void | zzfft (Int isign, Int n, Double scale, DComplex *x, DComplex *y, Double *table, Double *work, Int isys) |
| static void | scfft (Int isign, Int n, Float scale, Float *x, Complex *y, Float *table, Float *work, Int isys) |
scfft/dzfft computes the FFT of the real array x, and it stores the results in the complex array y. More... | |
| static void | scfft (Int isign, Int n, Double scale, Double *x, DComplex *y, Double *table, Double *work, Int isys) |
| static void | dzfft (Int isign, Int n, Double scale, Double *x, DComplex *y, Double *table, Double *work, Int isys) |
| static void | csfft (Int isign, Int n, Float scale, Complex *x, Float *y, Float *table, Float *work, Int isys) |
| static void | csfft (Int isign, Int n, Double scale, DComplex *x, Double *y, Double *table, Double *work, Int isys) |
| static void | zdfft (Int isign, Int n, Double scale, DComplex *x, Double *y, Double *table, Double *work, Int isys) |
| static void | ccfftm (Int isign, Int n, Int lot, Float scale, Complex *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys) |
ccfftm/zzfftm computes the FFT of each column of the complex matrix x, and stores the results in the columns of complex matrix y. More... | |
| static void | zzfftm (Int isign, Int n, Int lot, Double scale, DComplex *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys) |
| static void | scfftm (Int isign, Int n, Int lot, Float scale, Float *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys) |
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. More... | |
| static void | dzfftm (Int isign, Int n, Int lot, Double scale, Double *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys) |
| static void | csfftm (Int isign, Int n, Int lot, Float scale, Complex *x, Int ldx, Float *y, Int ldy, Float *table, Float *work, Int isys) |
| static void | zdfftm (Int isign, Int n, Int lot, Double scale, DComplex *x, Int ldx, Double *y, Int ldy, Double *table, Double *work, Int isys) |
| static void | ccfft2d (Int isign, Int n1, Int n2, Float scale, Complex *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys) |
| 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. More... | |
| static void | zzfft2d (Int isign, Int n1, Int n2, Double scale, DComplex *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys) |
| static void | scfft2d (Int isign, Int n1, Int n2, Float scale, Float *x, Int ldx, Complex *y, Int ldy, Float *table, Float *work, Int isys) |
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. More... | |
| static void | dzfft2d (Int isign, Int n1, Int n2, Double scale, Double *x, Int ldx, DComplex *y, Int ldy, Double *table, Double *work, Int isys) |
| static void | csfft2d (Int isign, Int n1, Int n2, Float scale, Complex *x, Int ldx, Float *y, Int ldy, Float *table, Float *work, Int isys) |
| static void | zdfft2d (Int isign, Int n1, Int n2, Double scale, DComplex *x, Int ldx, Double *y, Int ldy, Double *table, Double *work, Int isys) |
| static void | ccfft3d (Int isign, Int n1, Int n2, Int n3, Float scale, Complex *x, Int ldx, Int ldx2, Complex *y, Int ldy, Int ldy2, Float *table, Float *work, Int isys) |
| These routines compute the three-dimensional complex FFT of the complex matrix X, and store the results in the complex matrix Y. More... | |
| static void | zzfft3d (Int isign, Int n1, Int n2, Int n3, Double scale, DComplex *x, Int ldx, Int ldx2, DComplex *y, Int ldy, Int ldy2, Double *table, Double *work, Int isys) |
| static void | scfft3d (Int isign, Int n1, Int n2, Int n3, Float scale, Float *x, Int ldx, Int ldx2, Complex *y, Int ldy, Int ldy2, Float *table, Float *work, Int isys) |
| These are C++ wrapper functions for the 3D real-to-complex and complex-to-real transform routines in the SGI/Cray Scientific Library (SCSL). More... | |
| static void | dzfft3d (Int isign, Int n1, Int n2, Int n3, Double scale, Double *x, Int ldx, Int ldx2, DComplex *y, Int ldy, Int ldy2, Double *table, Double *work, Int isys) |
| static void | csfft3d (Int isign, Int n1, Int n2, Int n3, Float scale, Complex *x, Int ldx, Int ldx2, Float *y, Int ldy, Int ldy2, Float *table, Float *work, Int isys) |
| static void | zdfft3d (Int isign, Int n1, Int n2, Int n3, Double scale, DComplex *x, Int ldx, Int ldx2, Double *y, Int ldy, Int ldy2, Double *table, Double *work, Int isys) |
C++ Interface to the Sgi/Cray Scientific Library (SCSL)
These are C++ wrapper functions for the 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.
Warning:
Currently, the SCSL is available only on SGI machines;
|
static |
|
static |
These routines compute the Fast Fourier Transform (FFT) of the complex vector x, and store the result in vector y.
ccfft 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:
The output array is the FFT of the input array, using the following formula for the FFT:
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.
The table array stores the trigonometric tables used in calculation of the FFT. You must initialize table by calling the routine with isign = 0 prior to doing the transforms. If the value of the problem size, n, does not change, table does not have to be reinitialized.
In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:
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):
These examples use the table and workspace sizes appropriate to the Origin series.
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
ccfft: real. zzfft: double precision. Each element of the output array is multiplied by scale after taking the Fourier transform, as defined by the previous formula. Array of dimension (0:n-1). ccfft: complex array. zzfft: double complex array.
Input array of values to be transformed.
Output parameters:
ccfft: complex array. zzfft: double complex array. Output array of transformed values. The output array may be the same as the input array. In that case, the transform is done in place and the input array is overwritten with the transformed values. Real array; dimension 2*n+30.
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).
Real array; dimension 2*n.
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.
|
static |
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:
These routines compute the formula:
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).
These routines use a routine very much like ccfftm/zzfftm to do multiple FFTs first on all columns in an input matrix and then on all of the rows.
The table array stores factors of n1 and n2 and also trigonometric tables that are used in calculation of the FFT. This table must be initialized by calling the routine with isign = 0. If the values of the problem sizes, n1 and n2, do not change, the table does not have to be reinitialized.
In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:
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:
All examples here are for Origin series only.
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse transform as follows:
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.
Integer.
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).
Output parameters:
Real array; dimension (30+ 2 * n1) + (30 + 2 * n2).
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).
Real array; dimension 2 * (n1*n2).
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.
|
static |
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:
These routines compute the formula:
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).
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
Integer. Transform size in the first dimension. If n1 is not positive, the routine returns without computing a transform.
Integer. Transform size in the second dimension. If n2 is not positive, the routine returns without computing a transform.
Integer. Transform size in the third dimension. If n3 is not positive, the routine returns without computing a transform.
Scale factor. ccfft3d: real. zzfft3d: double precision.
Each element of the output array is multiplied by scale after taking the Fourier transform, as defined previously.
Array of dimension (0:ldx-1, 0:ldx2-1, 0:n3-1). ccfft3d: complex array. zzfft3d: double complex array.
Input array of values to be transformed.
Integer. The first dimension of x, as it was declared in the calling program (the leading dimension of x). ldx >= MAX(n1, 1).
Integer. The second dimension of x, as it was declared in the calling program. ldx2 >= MAX(n2, 1).
Integer. The first dimension of y, as it was declared in the calling program (the leading dimension of y). ldy >= MAX(n1, 1).
Integer. The second dimension of y, as it was declared in the calling program. ldy2 >= MAX(n2, 1).
Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.
isys = 0 or 1 depending on the amount of workspace the user can provide to the routine.
Output parameters:
Array of dimension (0:ldy-1, 0:ldy2-1, 0:n3-1). ccfft3d: complex array. zzfft3d: double complex array.
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.
Real array; dimension 30 + 2 * n1) + (30 + 2 * n2) + (30 + 2 * n3).
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).
Real array; dimension (n1 * n2 * n3).
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.
|
static |
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:
Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:
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.
These routines use decimation-in-frequency type FFT. It takes the FFT of the columns and vectorizes the operations along the rows of the matrix. Thus, the vector length in the calculations depends on the row size, and the strides for vector loads and stores are the leading dimensions, ldx and ldy.
The table array stores the trigonometric tables used in calculation of the FFT. You must initialize the table array by calling the routine with isign = 0 prior to doing the transforms. If the value of the problem size, n, does not change, table does not have to be reinitialized.
In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:
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:
Example 1: Initialize the TABLE array in preparation for doing an FFT of size 128. 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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
ccfftm: real. zzfftm: double precision. Each element of the output array is multiplied by scale factor after taking the Fourier transform, as defined previously. ccfftm: real array. zzfftm: double precision array. Input array of values to be transformed. Integer array of dimension (0:isys(0)). The first element of the array specifies how many more elements are in the array. Use isys to specify certain processor-specific parameters or options.
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 parameters:
Array of dimension (0:ldy-1, 0:lot-1). ccfftm: complex array. zzfftm: double complex array. Output array of transformed values. 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 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.
Real array; dimension (30 + 2n). 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).
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
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:
Then the output array is the FFT of the input array, using the following formula for the FFT:
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.
Notice in the preceding formula that there are n real input values, and n/2 + 1 complex output values. This property is characteristic of real-to-complex FFTs.
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:
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.
Consider the complex-to-real case. The effect of the computation is given by the preceding formula, but with X complex and Y real.
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:
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:
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.
The table array stores the trigonometric tables used in calculation of the FFT. This table must be initialized by calling the routine with isign = 0 prior to doing the transforms. The table does not have to be reinitialized if the value of the problem size, n, does not change. Because scfft and csfft use the same format for table, either can be used to initialize it (note that CCFFT uses a different table format).
In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared (assuming n > 0):
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:
These examples use the table and workspace sizes appropriate to Origin series.
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
scfft/dzfft returns without calculating the transform. scfft: real. dzfft: double precision. csfft: real. zdfft: double precision. Each element of the output array is multiplied by scale after taking the Fourier transform, as defined by the previous formula. scfft: real array of dimension (0:n-1). dzfft: double precision array of dimension (0:n-1). csfft: complex array of dimension (0:n/2). zdfft: double complex array of dimension (0:n/2). Integer array of dimension (0:isys(0)). Use isys to specify certain processor-specific parameters or options. The first element of the array specifies how many more elements are in the array.
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).
Output parameters:
Output array of transformed values. scfft: complex array of dimension (0:n/2). dzfft: double complex array of dimension (0:n/2). csfft: real array of dimension (0:n-1). zdfft: double precision array of dimension (0:n-1).
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.
Real array; dimension n+15.
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.
Real array; dimension n.
Work array used for intermediate calculations. Its address space must be different from that of the input and output arrays.
|
static |
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:
scfft2d computes the formula:
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:
Thus, you have to compute only (slightly more than) half of the output values, namely:
scfft2d uses a routine similar to scfftm to do a real-to-complex FFT on the columns, then uses a routine similar to ccfftm to do a complex-to-complex FFT on the rows.
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.
The table array stores factors of n1 and n2, and trigonometric tables that are used in calculation of the FFT. table must be initialized by calling the routine with isign = 0. table does not have to be reinitialized if the values of the problem sizes, n1 and n2, do not change.
In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared:
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:
The following examples are for Origin series only.
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
scfft2d returns without calculating a transform. scfft2d returns without calculating a transform. scfft2d: real. dzfft2d: double precision. csfft2d: real. zdfft2d: double precision. Each element of the output array is multiplied by scale factor after taking the Fourier transform, as defined previously. Array of dimension (0:ldx-1, 0:n2-1). scfft2d: real array. dzfft2d: double precision array. csfft2d: complex array. zdfft2d: double complex array.
Array of values to be transformed.
scfft2d, dzfft2d: ldx >= MAX(n1, 1). csfft2d, zdfft2d: ldx >= MAX(n1/2 + 1, 1). Integer.
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.
Integer array of dimension (0:isys(0)). The first element of the array specifies how many more elements are in the array. Use isys to specify certain processor-specific parameters or options.
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 parameters:
scfft2d: complex array. dzfft2d: double complex array. csfft2d: real array. zdfft2d: double precision array.
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.
Real array; dimension (15 + n1) + 2(15 + n2).
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).
Real array; dimension (n1 * n2).
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.
|
static |
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.
Warning:
Currently, the SCSL is available only on SGI machines;
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:
scfft3d computes the formula:
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
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:
Thus, you have to compute only (slightly more than) half out the output values, namely:
scfft3d uses a routine similar to scfftm to do multiple FFTs first on all columns of the input matrix, then uses a routine similar to ccfftm on all rows of the result, and then on all planes of that result. See scfftm and ccfftm for more information about the algorithms used.
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
Integer. Transform size in the first dimension. If n1 is not positive, scfft3d returns without computing a transform.
Integer. Transform size in the second dimension. If n2 is not positive, scfft3d returns without computing a transform.
Integer. Transform size in the third dimension. If n3 is not positive, scfft3d returns without computing a transform.
Scale factor. scfft3d: real. dzfft3d: double precision. csfft3d: real. zdfft3d: double precision. Each element of the output array is multiplied by scale after taking the Fourier transform, as defined previously.
Array of dimension (0:ldx-1, 0:ldx2-1, 0:n3-1). scfft3d: real array. dzfft3d: double precision array. csfft3d: complex array. zdfft3d: double complex array.
Array of values to be transformed.
Integer. The first dimension of x, as it was declared in the calling program (the leading dimension of x).
scfft3d, dzfft3d: ldx >= MAX(n1, 1). csfft3d, zdfft3d: ldx >= MAX(n1/2 + 1, 1).
Integer. The second dimension of x, as it was declared in the calling program. ldx2 >= MAX(n2, 1).
Integer. The first dimension of y, as it was declared in the calling program; that is, the leading dimension of y.
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.
Integer. The second dimension of y, as it was declared in the calling program. ldy2 >= MAX(n2, 1).
Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.
isys = 0 or 1 depending on the amount of workspace the user can provide to the routine.
Output parameters:
Array of dimension (0:ldy-1, 0:ldy2-1, 0:n3-1). scfft3d: complex array. dzfft3d: double complex array. csfft3d: real array. zdfft3d: double precision array.
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.
Real array; dimension (15 + n1) + 2(15 + n2) + 2(15 + n3).
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.
Real array; dimension n1 * n2 * n3.
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.
|
static |
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:
Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:
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.
Notice in the preceding formula that there are n real input values and (n/2) + 1 complex output values for each column. This property is characteristic of real-to-complex FFTs.
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:
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.
Consider the complex-to-real case. The effect of the computation is given by the preceding formula, but with X complex and Y real.
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:
Thus, in the complex-to-real routine, CSFFTM, the arrays can be dimensioned as follows:
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.
The table array contains the trigonometric tables used in calculation of the FFT. You must initialize this table by calling the routine with isign = 0 prior to doing the transforms. table does not have to be reinitialized if the value of the problem size, n, does not change.
In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared (for SCFFTM):
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:
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.
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.
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.
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.
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.
Input parameters:
Integer. Specifies whether to initialize the table array or to do the forward or inverse Fourier transform, as follows:
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.
scfftm or csfftm returns without computing a transforms. csfftm or scfftm returns without computing a transforms. scfftm: real. dzfftm: double precision. csfftm: real. zdfftm: double precision. Each element of the output array is multiplied by scale after taking the transform, as defined in the preceding formula. scfftm: real array. dzfftm: double precision array. csfftm: complex array. zdfftm: double complex array. scfftm, dzfftm: ldx >= MAX(n, 1). csfftm, zdfftm: ldx >= MAX(n/2 + 1, 1). scfftm, dzfftm: ldy >= MAX(n/2 + 1, 1). csfftm, zdfftm: ldy >= MAX(n, 1). Integer array of dimension (0:isys(0)). The first element of the array specifies how many more elements are in the array. Use isys to specify certain processor-specific parameters or options.
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 parameters:
Output array of transformed values. Dimension (0:ldy-1, 0:lot-1). scfftm: complex array. dzfftm: double complex array. csfftm: real array. zdfftm: double precision array.
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.
Real array; dimension (15 + n). Table of factors and trigonometric functions. This array must be initialized by a call to scfftm (or csfftm) with isign = 0.
If isign = 0, table is initialized to contain trigonometric tables needed to compute an FFT of length n.
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
|
static |
1.9.1