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

SCFT and DCFT--Complex Fourier Transform

These subroutines compute a set of m complex discrete n-point Fourier transforms of complex data.

Table 131. 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 transforms of the given sequences are 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, consisting of m sequences of length n. Specified as: an array of (at least) length 1+(n-1)inc1x+(m-1)inc2x, containing numbers of the data type indicated in Table 131.

inc1x

is the stride between the elements within each sequence in array X. Specified as: a fullword integer; inc1x > 0.

inc2x

is the stride between the first elements of the sequences in array X. (If m = 1, this argument is ignored.) Specified as: a fullword integer; inc2x > 0.

y

See On Return.

inc1y

is the stride between the elements within each sequence in array Y. Specified as: a fullword integer; inc1y > 0.

inc2y

is the stride between the first elements of each sequence in array Y. (If m = 1, this argument is ignored.) Specified as: a fullword integer; inc2y > 0.

n

is the length of each sequence to be transformed. Specified as: a fullword integer; n <= 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, as well as in the optionally-recoverable error 2030. For details, see Providing a Correct Transform Length to ESSL.

m

is the number of sequences to be transformed. Specified as: a fullword integer; m > 0.

isign

controls the direction of the transform, determining the sign Isign of the exponent of W_n, 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 131, 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 > 7 and naux1 >= (minimum value required for successful processing). To determine a sufficient value, use the processor-independent formulas. For values between 7 and 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, SCFT and DCFT 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, consisting of the results of the m discrete Fourier transforms, each of length n.

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

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 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. For optimal performance, the preferred value for inc1x and inc1y is 1. This implies that the sequences are stored with stride 1. The preferred value for inc2x and inc2y is n. This implies that sequences are stored one after another without any gap.

   It is possible to specify sequences in the transposed form--that is, as rows of a two-dimensional array. In this case, inc2x (or inc2y) = 1 and inc1x (or inc1y) is equal to the leading dimension of the array. One can specify either input, output, or both in the transposed form by specifying appropriate values for the stride parameters. For selecting optimal values of inc1x and inc1y for _CFT, you should use STRIDE--Determine the Stride Value for Optimal Performance in Specified Fourier Transform Subroutines. Example 1 in the STRIDE subroutine description explains how it is used for _CFT.

   If you specify the same array for X and Y, then inc1x and inc1y must be equal, and inc2x and inc2y must be equal. In this case, output overwrites input. If m = 1, the inc2x and inc2y values are not used by the subroutine. If you specify different arrays for X and Y, they must have no common elements; otherwise, results are unpredictable. See Concepts.

Processor-Independent Formulas for SCFT for NAUX1 and NAUX2

NAUX1 Formulas:

If n <= 8192, use naux1 = 20000.

If n > 8192, use naux1 = 20000+1.14n.

NAUX2 Formulas:

If n <= 8192, use naux2 = 20000.

If n > 8192, use naux2 = 20000+1.14n.

For the transposed case, where inc2x = 1 or inc2y = 1, and where n >= 252, add the following to the above storage requirements:

(n+256)(min(64, m)).

Processor-Independent Formulas for DCFT for NAUX1 and NAUX2

NAUX1 Formulas:

If n <= 2048, use naux1 = 20000.

If n > 2048, use naux1 = 20000+2.28n.

NAUX2 Formulas:

If n <= 2048, use naux2 = 20000.

If n > 2048, use naux2 = 20000+2.28n.

For the transposed case, where inc2x = 1 or inc2y = 1, and where n >= 252, add the following to the above storage requirements:

(2n+256)(min(64, m)).

Function

The set of m complex discrete n-point Fourier transforms of complex data in array X, with results going into array Y, is expressed as follows:

for:

k = 0, 1, ..., n-1

i = 1, 2, ..., m

where:

and where:

x_ji are elements of the sequences in array X.

y_ki are elements of the sequences in 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/n and isign being negative. See references [1], [3], [4], [19], 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 transforms.

Error Conditions

Resource Errors

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

Computational Errors

None

Input-Argument Errors

1. n > 37748736
2. inc1x, inc2x, inc1y, or inc2y <= 0
3. m <= 0
4. isign = 0
5. scale = 0.0
6. The subroutine has not been initialized with the present arguments.
7. The length of the transform in n is not an allowable value. Return code 1 is returned if error 2030 is recoverable.
8. naux1 <= 7
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 an input array X with a set of four short-precision complex sequences:

for j = 0, 1, ..., n-1 with n = 8, and the single frequencies k = 0, 1, 2, and 3. The arrays are declared as follows:

     COMPLEX*8  X(0:1023),Y(0:1023)
     REAL*8     AUX1(1693),AUX2(4096)

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 N   M ISIGN SCALE  AUX1   NAUX1  AUX2   NAUX2
           |    |   |     |   |   |     |   |   |   |     |     |       |     |       |
CALL SCFT(INIT, X , 1  ,  8 , Y , 1  ,  8 , 8 , 4 , 1 , SCALE, AUX1 , 1693 , AUX2 , 4096)

X contains the following four sequences:

(1.0000, 0.0000)   (1.0000,  0.0000)   (1.0000,  0.0000)   (1.0000,  0.0000)
(1.0000, 0.0000)   (0.7071,  0.7071)   (0.0000,  1.0000)  (-0.7071,  0.7071)
(1.0000, 0.0000)   (0.0000,  1.0000)  (-1.0000,  0.0000)   (0.0000, -1.0000)
(1.0000, 0.0000)  (-0.7071,  0.7071)   (0.0000, -1.0000)   (0.7071,  0.7071)
(1.0000, 0.0000)  (-1.0000,  0.0000)   (1.0000,  0.0000)  (-1.0000,  0.0000)
(1.0000, 0.0000)  (-0.7071, -0.7071)   (0.0000,  1.0000)   (0.7071, -0.7071)
(1.0000, 0.0000)   (0.0000, -1.0000)  (-1.0000,  0.0000)   (0.0000,  1.0000)
(1.0000, 0.0000)   (0.7071, -0.7071)   (0.0000, -1.0000)  (-0.7071, -0.7071)

Output

Y contains the following four sequences:

(8.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)
(0.0000, 0.0000)  (8.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)
(0.0000, 0.0000)  (0.0000, 0.0000)  (8.0000, 0.0000)  (0.0000, 0.0000)
(0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)  (8.0000, 0.0000)
(0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)
(0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)
(0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)
(0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)  (0.0000, 0.0000)

Example 2

This example shows an input array X with a set of four input spike sequences equal to the output of Example 1. This shows how you can compute the inverse of the transform in Example 1 by using a negative isign, giving as output the four sequences listed in the input for Example 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 INC1X INC2X Y INC1Y INC2Y N   M  ISIGN SCALE   AUX1   NAUX1  AUX2   NAUX2
           |    |   |     |   |   |     |   |   |    |     |      |       |     |       |
CALL SCFT(INIT, X , 1  ,  8 , Y , 1  ,  8 , 8 , 4 , -1 , SCALE , AUX1 , 1693 , AUX2 , 4096)

Output

Example 3

This example shows an input array X with a set of four short-precision complex sequences

for j = 0, 1, ..., n-1 with n = 12, and the single frequencies k = 0, 1, 2, and 3. Also, inc1x = inc1y = m and inc2x = inc2y = 1 to show how the input and output arrays can be stored in the transposed form. The arrays are declared as follows:

     COMPLEX*8  X (4,0:11),Y(4,0:11)
     REAL*8     AUX1(10000),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  N   M ISIGN SCALE  AUX1   NAUX1   AUX2   NAUX2
           |    |   |     |   |   |     |    |   |   |     |     |       |      |       |
CALL SCFT(INIT, X , 4  ,  1 , Y , 4  ,  1 , 12 , 4 , 1 , SCALE, AUX1 , 10000 , AUX2 , 10000)

X contains the following four sequences:

(1.0000, 0.0000)   (1.0000, 0.0000)   (1.0000, 0.0000)   (1.0000, 0.0000)
(1.0000, 0.0000)   (0.8660, 0.5000)   (0.5000, 0.8660)   (0.0000, 1.0000)
(1.0000, 0.0000)   (0.5000, 0.8660)  (-0.5000, 0.8660)  (-1.0000, 0.0000)
(1.0000, 0.0000)   (0.0000, 1.0000)  (-1.0000, 0.0000)  (0.0000, -1.0000)
(1.0000, 0.0000)  (-0.5000, 0.8660) (-0.5000, -0.8660)   (1.0000, 0.0000)
(1.0000, 0.0000)  (-0.8660, 0.5000)  (0.5000, -0.8660)   (0.0000, 1.0000)
(1.0000, 0.0000)  (-1.0000, 0.0000)   (1.0000, 0.0000)  (-1.0000, 0.0000)
(1.0000, 0.0000) (-0.8660, -0.5000)   (0.5000, 0.8660)  (0.0000, -1.0000)
(1.0000, 0.0000) (-0.5000, -0.8660)  (-0.5000, 0.8660)   (1.0000, 0.0000)
(1.0000, 0.0000)  (0.0000, -1.0000)  (-1.0000, 0.0000)   (0.0000, 1.0000)
(1.0000, 0.0000)  (0.5000, -0.8660) (-0.5000, -0.8660) (-1.0000,  0.0000)
(1.0000, 0.0000)  (0.8660, -0.5000)  (0.5000, -0.8660)  (0.0000, -1.0000)

Output

Y contains the following four sequences:

(12.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)  (12.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)  (12.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)  (12.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)
 (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)   (0.0000, 0.0000)

Example 4

This example shows an input array X with a set of four input spike sequences exactly equal to the output of Example 3. This shows how you can compute the inverse of the transform in Example 3 by using a negative isign, giving as output the four sequences listed in the input for Example 3. 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  N   M  ISIGN SCALE   AUX1  NAUX1  AUX2  NAUX2
           |    |   |     |   |   |     |    |   |    |     |      |      |     |      |
CALL SCFT(INIT, X , 4  ,  1 , Y , 4  ,  1 , 12 , 4 , -1 , SCALE , AUX1, 10000, AUX2, 10000)

Output

Example 5

This example shows how to compute a transform of a single long-precision complex sequence. It uses isign = 1 and scale = 1.0. The arrays are declared as follows:

     COMPLEX*16  X(0:7),Y(0:7)
     REAL*8      AUX1(26),AUX2(12)

The input in X is an impulse at zero, and the output in Y is constant for all frequencies. 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 N   M ISIGN SCALE   AUX1  NAUX1  AUX2  NAUX2
           |    |   |     |   |   |     |   |   |   |     |      |      |     |      |
CALL DCFT(INIT, X , 1  ,  0 , Y , 1  ,  0 , 8 , 1 , 1 , SCALE , AUX1 , 26  , AUX2 , 12)

X contains the following sequence:

(1.0000, 0.0000)
(0.0000, 0.0000)
(0.0000, 0.0000)
(0.0000, 0.0000)
(0.0000, 0.0000)
(0.0000, 0.0000)
(0.0000, 0.0000)
(0.0000, 0.0000)

Output

(1.0000, 0.0000)
(1.0000, 0.0000)
(1.0000, 0.0000)
(1.0000, 0.0000)
(1.0000, 0.0000)
(1.0000, 0.0000)
(1.0000, 0.0000)
(1.0000, 0.0000)