FAST_2DFT

Computes the Discrete Fourier Transform (2DFT) of a rank-2 complex array, x.

Required Arguments

No required arguments; pairs of optional arguments are required. These pairs are forward_in and forward_out or inverse_in and inverse_out.

Optional Arguments

forward_in = x (Input)
Stores the input complex array of rank-2 to be transformed.

forward_out = y (Output)
Stores the output complex array of rank-2 resulting from the transform.

inverse_in = y (Input)
Stores the input complex array of rank-2 to be inverted.

inverse_out = x (Output)
Stores the output complex array of rank-2 resulting from the inverse transform.

mdata = m (Input)
Uses the sub-array in first dimension of size m for the numbers.
Default value: m = size(x,1).

ndata = n (Input)
Uses the sub-array in the second dimension of size n for the numbers.
Default value: n = size(x,2).

ido = ido (Input/Output)
Integer flag that directs user action. Normally, this argument is used only when the working variables required for the transform and its inverse are saved in the calling program unit. Computing the working variables and saving them in internal arrays within fast_2dft is the default. This initialization step is expensive.

There is a two-step process to compute the working variables just once. Example 3 illustrates this usage. The general algorithm for this usage is to enter fast_2dft with
ido = 0. A return occurs thereafter with ido < 0. The optional rank-1 complex array w(:) with size(w) >= -ido must be re-allocated. Then, re-enter fast_2dft. The next return from fast_2dft has the output value ido = 1. The variables required for the transform and its inverse are saved in w(:). Thereafter, when the routine is entered with ido = 1 and for the same values of m and n, the contents of w(:) will be used for the working variables. The expensive initialization step is avoided. The optional arguments “ido=” and “work_array=” must be used together.

work_array = w(:) (Output/Input)
Complex array of rank-1 used to store working variables and values between calls to fast_2dft. The value for size(w) must be at least as large as the value - ido for the value of ido < 0.

iopt = iopt(:) (Input/Output)
Derived type array with the same precision as the input array; used for passing optional data to fast_2dft. The options are as follows:

Packaged Options for FAST_2DFT
Option Prefix = ?	Option Name	Option Value
c_, z_	fast_2dft_scan_for_NaN	1
c_, z_	fast_2dft_near_power_of_2	2
c_, z_	fast_2dft_scale_forward	3
c_, z_	fast_2dft_scale_inverse	4

iopt(IO) = ?_options(?_fast_2dft_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that
isNaN(x(i,j)) ==.true.
See the isNaN() function, Chapter 10 .
Default: Does not scan for NaNs.

iopt(IO) = ?_options(?_fast_2dft_near_power_of_2, ?_dummy)
Nearest powers of 2 ≥ m and ≥ n are returned as an outputs in iopt(IO + 1)%idummy and iopt(IO + 2)%idummy.

iopt(IO) = ?_options(?_fast_2dft_scale_forward, real_part_of_scale)

iopt(IO+1) = ?_options(?_dummy, imaginary_part_of_scale)
Complex number defined by the factor
cmplx(real_part_of_scale, imaginary_part_of_scale) is
multiplied by the forward transformed array.
Default value is 1.

iopt(IO) = ?_options(?_fast_2dft_scale_inverse, real_part_of_scale)

FORTRAN 90 Interface

Generic: None

Specific: The specific interface names are S_FAST_2DFT, D_FAST_2DFT, C_FAST_2DFT, and Z_FAST_2DFT.

Description

The fast_2dft routine is a Fortran 90 version of the FFT suite of IMSL (1994, pp. 772-776).

Fatal and Terminal Messages

See the messages.gls file for error messages for FAST_2DFT. These error messages are numbered 670 -680; 720 -730.

Example 1: Transforming an Array of Random Complex Numbers

An array of random complex numbers is obtained. The transform of the numbers is inverted and the final results are compared with the input array.

use fast_2dft_int

use rand_int

implicit none

! This is Example 1 for FAST_2DFT.

integer, parameter :: n=24

integer, parameter :: m=40

real(kind(1e0)) :: err, one=1e0

complex(kind(1e0)), dimension(n,m) :: a, b, c

! Generate a random complex sequence.

a=rand(a); c=a

! Transform and then invert the transform.

call c_fast_2dft(forward_in=a, &

forward_out=b)

call c_fast_2dft(inverse_in=b, &

inverse_out=a)

! Check that inverse(transform(sequence)) = sequence.

err = maxval(abs(c-a))/maxval(abs(c))

if (err <= sqrt(epsilon(one))) then

write (*,*) 'Example 1 for FAST_2DFT is correct.'

end if

end

Output

Example 1 for FAST_2DFT is correct.

Additional Examples

.p>.F90CH3.DOC!EX2_FAST_2DFT;Example 2: Cyclical 2D Data with a Linear Trend

This set of data is sampled from a function x(s, t) = a + bs + ct + y(s, t) , where y(s, t) is an harmonic series. The independent variables are normalized as -1 ≤ s ≤ 1 and -1 ≤ t ≤ 1. Thus, the data is said to have cyclical components plus a linear trend. As a first step, the linear terms are effectively removed from the data using the least-squares system solver . Then, the residuals are transformed and the resulting frequencies are analyzed.

use fast_2dft_int

use lin_sol_lsq_int

use sort_real_int

use rand_int

implicit none

! This is Example 2 for FAST_2DFT.

integer i

integer, parameter :: n=8, k=15

integer ip(n*n), order(k)

real(kind(1e0)), parameter :: one=1e0, two=2e0, zero=0e0

real(kind(1e0)) delta_t

real(kind(1e0)) rn(3), s(n), t(n), temp(n*n), new_order(k)

complex(kind(1e0)) a, b, c, a_trend(n*n,3), b_trend(n*n,1), &

f(n,n), r(n,n), x(n,n), x_trend(3,1)

complex(kind(1e0)), dimension(n,n) :: g=zero, h=zero

! Generate random data for planar trend.

rn = rand(rn)

a = rn(1)

b = rn(2)

c = rn(3)

! Generate the frequency components of the harmonic series.

! Non-zero random amplitudes given on two edges of the square domain.

g(1:,1)=rand(g(1:,1))

g(1,1:)=rand(g(1,1:))

! Invert 'g' into the harmonic series 'h' in time domain.

call c_fast_2dft(inverse_in=g, inverse_out=h)

! Compute sampling interval.

delta_t = two/n

s = (/(-one + (i-1)*delta_t, i=1,n)/)

t = (/(-one + (i-1)*delta_t, i=1,n)/)

! Make up data set as a linear trend plus harmonics.

x = a + b*spread(s,dim=2,ncopies=n) + &

c*spread(t,dim=1,ncopies=n) + h

! Define least-squares matrix data for a planar trend.

a_trend(1:,1) = one

a_trend(1:,2) = reshape(spread(s,dim=2,ncopies=n),(/n*n/))

a_trend(1:,3) = reshape(spread(t,dim=1,ncopies=n),(/n*n/))

b_trend(1:,1) = reshape(x,(/n*n/))

! Solve for a linear trend.

call lin_sol_lsq(a_trend, b_trend, x_trend)

! Compute harmonic residuals.

r = x - reshape(matmul(a_trend,x_trend),(/n,n/))

! Transform harmonic residuals.

call c_fast_2dft(forward_in=r, forward_out=f)

ip = (/(i,i=1,n**2)/)

! Sort the magnitude of the transform.

call s_sort_real(-(abs(reshape(f,(/n*n/)))), &

temp, iperm=ip)

! The dominant frequencies are output in ip(1:k).

! Sort these values to compare with the original frequency order.

call s_sort_real(real(ip(1:k)), new_order)

order(1:n) = (/(i,i=1,n)/)

order(n+1:k) = (/((i-n)*n+1,i=n+1,k)/)

! Check the results.

if (count(order /= int(new_order)) == 0) then

write (*,*) 'Example 2 for FAST_2DFT is correct.'

end if

end

Output

Example 2 for FAST_2DFT is correct.

Example 3: Several 2D Transforms with Initialization

In this example, the optional arguments ido and work_array are used to save working variables in the calling program unit. This results in maximum efficiency of the transform and its inverse since the working variables do not have to be precomputed following each entry to routine fast_2dft.

use fast_2dft_int

implicit none

! This is Example 3 for FAST_2DFT.

integer i, j

integer, parameter :: n=256

real(kind(1e0)), parameter :: one=1e0, zero=0e0

real(kind(1e0)) r(n,n), err

complex(kind(1e0)) a(n,n), b(n,n), c(n,n)

! The value of the array size for work(:) is computed in the

! routine fast_dft as a first step.

integer ido_value

complex(kind(1e0)), allocatable :: work(:)

! Fill in value one for points inside the circle with r=64.

a = zero

r = reshape((/(((i-n/2)**2 + (j-n/2)**2, i=1,n), &

j=1,n)/),(/n,n/))

where (r <= (n/4)**2) a = one

c = a

! Transform and then invert the sequence using the pre-computed

! working values.

ido_value = 0

if(allocated(work)) deallocate(work)

! Allocate the space required for work(:).

if (ido_value <= 0) allocate(work(-ido_value))

! Transform the image and then invert it back.

call c_fast_2dft(forward_in=a, &

forward_out=b, IDO=ido_value, work_array=work)

if (ido_value == 1) exit

end do

call c_fast_2dft(inverse_in=b, &

inverse_out=a, IDO=ido_value, work_array=work)

! Deallocate the space used for work(:).

if (allocated(work)) deallocate(work)

! Check that inverse(transform(image)) = image.

err = maxval(abs(c-a))/maxval(abs(c))

if (err <= sqrt(epsilon(one))) then

write (*,*) 'Example 3 for FAST_2DFT is correct.'

end if

end

Output

Example 3 for FAST_2DFT is correct.

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260