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.
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.
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.
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.
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).
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).
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.
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:
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.
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 outputs in iopt(IO + 1)%idummy and iopt(IO + 2)%idummy.
Nearest powers of 2≥m and ≥n are returned as 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.
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)
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 inverse transformed array.
Default value is 1.
Complex number defined by the factor cmplx(real_part_of_scale, imaginary_part_of_scale) is multiplied by the inverse transformed array.
Default value is 1.
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.
Examples
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.
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
do
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.
