IBM Books

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

SCRFT2 and DCRFT2--Complex-to-Real Fourier Transform in Two Dimensions

These subroutines compute the two-dimensional discrete Fourier transform of complex conjugate even data in a two-dimensional array.

Table 138. Data Types

X Y, scale Subroutine
Short-precision complex Short-precision real SCRFT2
Long-precision complex Long-precision real DCRFT2
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

Fortran CALL SCRFT2 (init, x, inc2x, y, inc2y, n1, n2, isign, scale, aux1, naux1, aux2, naux2, aux3, naux3)

CALL DCRFT2 (init, x, inc2x, y, inc2y, n1, n2, isign, scale, aux1, naux1, aux2, naux2)

C and C++ scrft2 (init, x, inc2x, y, inc2y, n1, n2, isign, scale, aux1, naux1, aux2, naux2, aux3, naux3);

dcrft2 (init, x, inc2x, y, inc2y, n1, n2, isign, scale, aux1, naux1, aux2, naux2);

PL/I CALL SCRFT2 (init, x, inc2x, y, inc2y, n1, n2, isign, scale, aux1, naux1, aux2, naux2, aux3, naux3);

CALL DCRFT2 (init, x, inc2x, y, inc2y, n1, n2, isign, scale, aux1, naux1, aux2, naux2);

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 n2 columns of data to be transformed. Due to complex conjugate symmetry, the input consists of only the first ((n1)/2)+1 rows of the array--that is, xj1,j2, j1 = 0, 1, ..., (n1)/2, j2 = 0, 1, ..., n2-1. The data in each column is stored with stride 1.

Specified as: an inc2x by (at least) n2 array, containing numbers of the data type indicated in Table 138. This array must be aligned on a doubleword boundary.

inc2x
is the leading dimension (stride between columns) of array X. Specified as: a fullword integer; inc2x >= ((n1)/2)+1.

y
See On Return.

inc2y
is the leading dimension (stride between the columns) of array Y. Specified as: a fullword integer; inc2y >= n1+2.

n1
is the number of rows of data--that is, the length of the columns in array Y involved in the computation. The length of the columns in array X are (n1)/2+1. 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 number of columns of data--that is, the length of the rows in arrays X and Y involved in the computation. 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 Wn1 and Wn2, 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 138, 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 aux. 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, SCRFT2 and DCRFT2 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.

aux3
this argument is provided for migration purposes only and is ignored.

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

naux3
this argument is provided for migration purposes only and is ignored.

Specified as: a fullword integer.

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 the array Y, containing n1 rows and n2 columns of results of the real discrete Fourier transform of X. The data in each column of Y is stored with stride 1.

Returned as: an inc2y by (at least) n2 array, containing numbers of the data type indicated in Table 138. 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 (2)(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. For selecting optimal strides (or leading dimensions inc2x and inc2y) for your input and output arrays, you should use STRIDE--Determine the Stride Value for Optimal Performance in Specified Fourier Transform Subroutines. Example 6 in the STRIDE subroutine description explains how it is used for these subroutines.
  5. Be sure to align array Y on a doubleword boundary, and specify an even number for inc2y, if possible.

Processor-Independent Formulas for SCRFT2 for NAUX1 and NAUX2

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

NAUX1 Formulas
If max(n1/2, n2) <= 8192, use naux1 = 45000.
If max(n1/2, n2) > 8192, use naux1 = 40000+0.82n1+1.14n2.

NAUX2 Formulas
If n1 <= 16384 and n2 < 252, use naux2 = 20000.
If n1 > 16384 and n2 < 252, use naux2 = 20000+0.57n1.
If n2 >= 252, add the following to the above storage requirements:
(n2+256)(1.14+s)
where s = min(64, 1+n1/2).

Processor-Independent Formulas for DCRFT2 for NAUX1 and NAUX2

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

NAUX1 Formulas
If n <= 2048, use naux1 = 42000.
If n > 2048, use naux1 = 40000+1.64n1+2.28n2,
where n = max(n1/2, n2).

NAUX2 Formulas
If n1 <= 4096 and n2 < 252, use naux2 = 20000.
If n1 > 4096 and n2 < 252, use naux2 = 20000+1.14n1.
If n2 >= 252, add the following to the above storage requirements:
((2)n2+256) (2.28+s)
where s = min(64, 1+n1/2).

Function

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



Two-Dimensional FFT Graphic

for:

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

where:



Two-Dimensional FFT Graphic

and where:

xj1,j2 are elements of array X.
yk1,k2 are elements of array Y.
Isign is + or - (determined by argument isign).
scale is a scalar value.

Because of the complex conjugate symmetry, the output in array Y is real. 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. inc2x < (n1)/2+1
  4. inc2y < n1+2
  5. scale = 0.0
  6. isign = 0
  7. The subroutine has not been initialized with the present arguments.
  8. 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.
  9. naux1 is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.
  10. 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. The arrays are declared as follows: 

     REAL*4     Y(0:13,0:7)
     COMPLEX*8  X(0:6,0:7)
     REAL*8     AUX1(1000), AUX2(1000), AUX3(1)

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 INC2X Y  INC2Y N1   N2 ISIGN SCALE   AUX1  NAUX1   AUX2  NAUX2   AUX3  NAUX3
             |    |   |   |    |   |    |    |     |      |      |      |      |      |      |
CALL SCRFT2(INIT, X , 7 , Y , 14 , 12 , 8 , -1 , SCALE , AUX1 , 1000 , AUX2 , 1000 , AUX3 ,  0 )

INIT = 1(for initialization)
INIT = 0(for computation)
SCALE = 1.0/96.0
X is an array with 7 rows and 8 columns with (1.0, 0.0) in all locations.

Output
        *                                        *
        | 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
Y    =  | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 |
        |  .    .    .    .    .    .    .    .  |
        |  .    .    .    .    .    .    .    .  |
        *                                        *

Example 2

This example shows another transform computation with different data using the same initialized array AUX1 in Example 1.

Call Statement and Input


            INIT X INC2X Y  INC2Y N1   N2 ISIGN SCALE   AUX1   NAUX1  AUX2   NAUX2  AUX3  NAUX3
             |   |   |   |    |   |    |    |     |      |       |     |       |     |      |
CALL SCRFT2( 0 , X , 7 , Y , 14 , 12 , 8 , -1 , SCALE , AUX1 , 1000 , AUX2 , 1000 , AUX3 ,  0 )

SCALE    =  1.0/96.0

X is an array with 7 rows and 8 columns having (96.0, 0.0) in location X(0,0) and (0.0, 0.0) in all others.

Output
        *                                        *
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
Y    =  | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        |  .    .    .    .    .    .    .    .  |
        |  .    .    .    .    .    .    .    .  |
        *                                        *

Example 3

This example shows the same array being used for input and output. The arrays are declared as follows: 

     REAL*8      Y(0:17,0:7)
     COMPLEX*16  X(0:8,0:7)
     REAL*8      AUX1(1000), AUX2(1000), AUX3(1)

Arrays X and Y are made equivalent by the following statement, making them occupy the same storage. 

     EQUIVALENCE (X,Y)

This requires inc2y = 2(inc2x). 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 INC2X Y  INC2Y N1  N2 ISIGN SCALE   AUX1   NAUX1  AUX2   NAUX2  AUX3  NAUX3
             |    |   |   |    |    |   |   |     |      |       |      |      |     |      |
CALL DCRFT2(INIT, X , 9 , Y , 18 , 16 , 8 , 1 , SCALE , AUX1 , 1000 , AUX2 , 1000 , AUX3 ,  0 )

INIT = 1(for initialization)
INIT = 0(for computation)
SCALE = 1.0

X is an array with 9 rows and 8 columns having (1.0, 1.0) in location X(4,2) and (0.0, 0.0) in all others.

Output
        *                                                *
        |  2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0 |
        |  2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0 |
        | -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0 |
        | -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0 |
        |  2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0 |
        |  2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0 |
        | -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0 |
        | -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0 |
Y    =  |  2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0 |
        |  2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0 |
        | -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0 |
        | -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0 |
        |  2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0 |
        |  2.0  -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0 |
        | -2.0  -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0 |
        | -2.0   2.0   2.0  -2.0  -2.0   2.0   2.0  -2.0 |
        |   .     .     .     .     .     .     .     .  |
        |   .     .     .     .     .     .     .     .  |
        *                                                *

[ Top of Page | Previous Page | Next Page | Table of Contents | Index ]