Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

SCFT2 and DCFT2--Complex Fourier Transform in Two Dimensions

These subroutines compute the two-dimensional discrete Fourier transform of complex data.

Table 136. Data Types

Note:

Two invocations of this subroutine are necessary: one to prepare the working storage for the subroutine, and the other to perform the computations.

Syntax

On Entry

init

is a flag, where:

If init <> 0, trigonometric functions and other parameters, depending on arguments other than x, are computed and saved in aux1. The contents of x and y are not used or changed.

If init = 0, the discrete Fourier transform of the given array is computed. The only arguments that may change after initialization are x, y, and aux2. All scalar arguments must be the same as when the subroutine was called for initialization with init <> 0.

Specified as: a fullword integer. It can have any value.

x

is the array X, containing the two-dimensional data to be transformed, where each element x_j1,j2, using zero-based indexing, is stored in X(j1(inc1x)+j2(inc2x)) for j1 = 0, 1, ..., n1-1 and j2 = 0, 1, ..., n2-1.

Specified as: an array of (at least) length 1+(n1-1)inc1x+(n2-1)inc2x, containing numbers of the data type indicated in Table 136. This array must be aligned on a doubleword boundary, and:

If inc1x = 1, the input array is stored in normal form, and inc2x >= n1.

If inc2x = 1, the input array is stored in transposed form, and inc1x >= n2.

See Notes for more details.

inc1x

is the stride between the elements in array X for the first dimension.

If the array is stored in the normal form, inc1x = 1.

If the array is stored in the transposed form, inc1x is the leading dimension of the array and inc1x >= n2.

Specified as: a fullword integer; inc1x > 0. If inc2x = 1, then inc1x >= n2.

inc2x

is the stride between the elements in array X for the second dimension.

If the array is stored in the transposed form, inc2x = 1.

If the array is stored in the normal form, inc2x is the leading dimension of the array and inc2x >= n1.

Specified as: a fullword integer; inc2x > 0. If inc1x = 1, then inc2x >= n1.

y

See On Return.

inc1y

is the stride between the elements in array Y for the first dimension.

If the array is stored in the normal form, inc1y = 1.

If the array is stored in the transposed form, inc1y is the leading dimension of the array and inc1y >= n2.

Specified as: a fullword integer; inc1y > 0. If inc2y = 1, then inc1y >= n2.

inc2y

is the stride between the elements in array Y for the second dimension.

If the array is stored in the transposed form, inc2y = 1.

If the array is stored in the normal form, inc2y is the leading dimension of the array and inc2y >= n1.

Specified as: a fullword integer; inc2y > 0. If inc1y = 1, then inc2y >= n1.

n1

is the length of the first dimension of the two-dimensional data in the array to be transformed. Specified as: a fullword integer; n1 <= 37748736 and must be one of the values listed in Acceptable Lengths for the Transforms. For all other values specified less than 37748736, you have the option of having the next larger acceptable value returned in this argument. For details, see Providing a Correct Transform Length to ESSL.

n2

is the length of the second dimension of the two-dimensional data in the array to be transformed. Specified as: a fullword integer; n2 <= 37748736 and must be one of the values listed in Acceptable Lengths for the Transforms. For all other values specified less than 37748736, you have the option of having the next larger acceptable value returned in this argument. For details, see Providing a Correct Transform Length to ESSL.

isign

controls the direction of the transform, determining the sign Isign of the exponents of W_n1 and W_n2, where:

If isign = positive value, Isign = + (transforming time to frequency).

If isign = negative value, Isign = - (transforming frequency to time).

Specified as: a fullword integer; isign > 0 or isign < 0.

scale

is the scaling constant scale. See Function for its usage. Specified as: a number of the data type indicated in Table 136, where scale > 0.0 or scale < 0.0.

aux1

is the working storage for this subroutine, where:

If init <> 0, the working storage is computed.

If init = 0, the working storage is used in the computation of the Fourier transforms.

Specified as: an area of storage, containing naux1 long-precision real numbers.

naux1

is the number of doublewords in the working storage specified in aux1. Specified as: a fullword integer; naux1 >= (minimum value required for successful processing). To determine a sufficient value, use the processor-independent formulas. For all other values specified less than the minimum value, you have the option of having the minimum value returned in this argument. For details, see Using Auxiliary Storage in ESSL.

aux2

has the following meaning:

If naux2 = 0 and error 2015 is unrecoverable, aux2 is ignored.

Otherwise, it is the working storage used by this subroutine, which is available for use by the calling program between calls to this subroutine.

Specified as: an area of storage, containing naux2 long-precision real numbers. On output, the contents are overwritten.

naux2

is the number of doublewords in the working storage specified in aux2. Specified as: a fullword integer, where:

If naux2 = 0 and error 2015 is unrecoverable, SCFT2 and DCFT2 dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux2 >= (minimum value required for successful processing). To determine a sufficient value, use the processor-independent formulas. For all other values specified less than the minimum value, you have the option of having the minimum value returned in this argument. For details, see Using Auxiliary Storage in ESSL.

On Return

y

has the following meaning, where:

If init <> 0, this argument is not used, and its contents remain unchanged.

If init = 0, this is array Y, containing the elements resulting from the two-dimensional discrete Fourier transform of the data in X. Each element y_k1,k2, using zero-based indexing, is stored in Y(k1(inc1y)+k2(inc2y)) for k1 = 0, 1, ..., n1-1 and k2 = 0, 1, ..., n2-1.

Returned as: an array of (at least) length 1+(n1-1)inc1y+(n2-1)inc2y, containing numbers of the data type indicated in Table 136. This array must be aligned on a doubleword boundary, and:

If inc1y = 1, the output array is stored in normal form, and inc2y >= n1.

If inc2y = 1, the output array is stored in transposed form, and inc1y >= n2.

See Notes for more details.

aux1

is the working storage for this subroutine, where:

If init <> 0, it contains information ready to be passed in a subsequent invocation of this subroutine.

If init = 0, its contents are unchanged.

Returned as: the contents are not relevant.

Notes

1. aux1 should not be used by the calling program between program calls to this subroutine with init <> 0 and init = 0. However, it can be reused after intervening calls to this subroutine with different arguments.
2. When using the ESSL SMP library, for optimal performance, the number of threads specified should be the same for init <> 0 and init = 0.
3. If you specify the same array for X and Y, then inc1x must equal inc1y, and inc2x must equal inc2y. In this case, output overwrites input. If you specify different arrays X and Y, they must have no common elements; otherwise, results are unpredictable. See Concepts.
4. By appropriately specifying the inc arguments, this subroutine allows you to specify that it should use one of two forms of its arrays, the normal untransposed form or the transposed form. As a result, you do not have to move any data. Instead, the subroutine performs the adjustments for you. Also, either the input array or the output array can be in transposed form. The FFT computation is symmetrical with respect to n1 and n2. They can be interchanged without the loss of generality. If they are interchanged, an array that is stored in the normal form appears as an array stored in the transposed form and vise versa. If, for performance reasons, the forms of the input and output arrays are different, then the input array should be specified in the normal form, and the output array should be specified in the transposed form. This can always be done by interchanging n1 and n2.
5. Although the inc arguments for each array can be arbitrary, in most cases, one of the inc arguments is 1 for each array. If inc1 = 1, the array is stored in normal form; that is, the first dimension of the array is along the columns. In this case, inc2 is the leading dimension of the array and must be at least n1. Conversely, if inc2 = 1, the array is stored in the transposed form; that is, the first dimension of the array is along the rows. In this case, inc1 is the leading dimension of the array and must be at least n2. The rows of the arrays are accessed with a stride that equals the leading dimension of the array. To minimize cache interference in accessing a row, an optimal value should be used for the leading dimension of the array. You should use STRIDE--Determine the Stride Value for Optimal Performance in Specified Fourier Transform Subroutines to determine this optimal value. Example 4 in the STRIDE subroutine description explains how it is used to find either inc1 or inc2.

Processor-Independent Formulas for SCFT2 for NAUX1 and NAUX2

The required values of naux1 and naux2 depend on n1 and n2.

NAUX1 Formulas

If max(n1, n2) <= 8192, use naux1 = 40000.

If max(n1, n2) > 8192, use naux1 = 40000+1.14(n1+n2).

NAUX2 Formulas

If max(n1, n2) < 252, use naux2 = 20000.

If max(n1, n2) >= 252, use naux2 = 20000+(r+256)(s+1.14), where r = max(n1, n2) and s = min(64, n1, n2).

Processor-Independent Formulas for DCFT2 for NAUX1 and NAUX2

The required values of naux1 and naux2 depend on n1 and n2.

NAUX1 Formulas

If max(n1, n2) <= 2048, use naux1 = 40000.

If max(n1, n2) > 2048, use naux1 = 40000+2.28(n1+n2).

NAUX2 Formulas

If max(n1, n2) < 252, use naux2 = 20000.

If max(n1, n2) >= 252, use naux2 = 20000+(2r+256)(s+2.28), where r = max(n1, n2) and s = min(64, n1, n2).

Function

The two-dimensional discrete Fourier transform of complex data in array X, with results going into array Y, is expressed as follows:

for:

k1 = 0, 1, ..., n1-1

k2 = 0, 1, ..., n2-1

where:

and where:

x_j1,j2 are elements of array X.

y_k1,k2 are elements of array Y.

Isign is + or - (determined by argument isign).

scale is a scalar value.

For scale = 1.0 and isign being positive, you obtain the discrete Fourier transform, a function of frequency. The inverse Fourier transform is obtained with scale = 1.0/((n1)(n2)) and isign being negative. See references [1], [4], and [20].

Two invocations of this subroutine are necessary:

1. With init <> 0, the subroutine tests and initializes arguments of the program, setting up the aux1 working storage.
2. With init = 0, the subroutine checks that the initialization arguments in the aux1 working storage correspond to the present arguments, and if so, performs the calculation of the Fourier transform.

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux2 = 0, and unable to allocate work area.

Computational Errors

None

Input-Argument Errors

1. n1 > 37748736
2. n2 > 37748736
3. inc1x|inc2x|inc1y|inc2y <= 0
4. scale = 0.0
5. isign = 0
6. The subroutine has not been initialized with the present arguments.
7. The length of one of the transforms in n1 or n2 is not an allowable value. Return code 1 is returned if error 2030 is recoverable.
8. naux1 is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.
9. Error 2015 is recoverable or naux2<>0, and naux2 is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example 1

This example shows how to compute a two-dimensional transform where both input and output are stored in normal form (inc1x = inc1y = 1). Also, inc2x = inc2y so the same array can be used for both input and output. The arrays are declared as follows:

     COMPLEX*8  X(6,8),Y(6,8)
     REAL*8     AUX1(20000), AUX2(10000)

Arrays X and Y are made equivalent by the following statement, making them occupy the same storage: EQUIVALENCE (X,Y). First, initialize AUX1 using the calling sequence shown below with INIT <> 0. Then use the same calling sequence with INIT = 0 to do the calculation.

Call Statement and Input

           INIT  X INC1X INC2X Y INC1Y INC2Y N1  N2 ISIGN SCALE  AUX1  NAUX1   AUX2  NAUX2
            |    |   |     |   |   |     |   |   |    |     |      |     |      |      |
CALL SCFT2(INIT, X , 1  ,  6 , Y , 1  ,  6 , 6 , 8 ,  1 , SCALE, AUX1, 20000 , AUX2, 10000)

Output

Y is an array with 6 rows and 8 columns having (48.0, 0.0) in location Y(1,1) and (0.0, 0.0) in all others.

Example 2

This example shows how to compute a two-dimensional inverse Fourier transform. For this example, X is stored in normal untransposed form (inc1x = 1), and Y is stored in transposed form (inc2y = 1). The arrays are declared as follows:

     COMPLEX*16  X(6,8),Y(8,6)
     REAL*8      AUX1(20000), AUX2(10000)

First, initialize AUX1 using the calling sequence shown below with INIT <> 0. Then use the same calling sequence with INIT = 0 to do the calculation.

Call Statement and Input

           INIT  X INC1X INC2X Y INC1Y INC2Y N1  N2 ISIGN SCALE  AUX1   NAUX1   AUX2   NAUX2
            |    |   |     |   |   |     |   |   |    |     |     |       |      |       |
CALL DCFT2(INIT, X , 1  ,  6 , Y , 8  ,  1 , 6 , 8 , -1 , SCALE, AUX1 , 20000 , AUX2 , 10000)

Output

Y is an array with 8 rows and 6 columns with (1.0, 0.0) in all locations.