FAST_DFT
Computes the Discrete Fourier Transform (DFT) of a rank-1 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-1 to be transformed.
Stores the input complex array of rank-1 to be transformed.
forward_out = y (Output)
Stores the output complex array of rank-1 resulting from the transform.
Stores the output complex array of rank-1 resulting from the transform.
inverse_in = y (Input)
Stores the input complex array of rank-1 to be inverted.
Stores the input complex array of rank-1 to be inverted.
inverse_out = x (Output)
Stores the output complex array of rank-1 resulting from the inverse transform.
Stores the output complex array of rank-1 resulting from the inverse transform.
ndata = n (Input)
Uses the sub-array of size n for the numbers.
Default value: n = size(x).
Uses the sub-array of size n for the numbers.
Default value: n = size(x).
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_dft 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_dft 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_dft 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_dft. The next return from fast_dft 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 value of 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_dft. The value for size(w) must be at least as large as the value – ido for the value of ido < 0.
Complex array of rank-1 used to store working variables and values between calls to fast_dft. 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_dft. The options are as follows:
Derived type array with the same precision as the input array; used for passing optional data to fast_dft. The options are as follows:
Packaged Options for FAST_DFT | ||
|---|---|---|
Option Prefix = ? | Option Name | Option Value |
c_, z_ | fast_dft_scan_for_NaN | 1 |
c_, z_ | fast_dft_near_power_of_2 | 2 |
c_, z_ | fast_dft_scale_forward | 3 |
c_, z_ | Fast_dft_scale_inverse | 4 |
iopt(IO) = ?_options(?_fast_dft_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that isNaN(x(i)) ==.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)) ==.true..
See the isNaN() function, Chapter 10.
Default: Does not scan for NaNs.
iopt(IO) = ?_options(?_fast_dft_near_power_of_2, ?_dummy)
Nearest power of 2 ≥ n is returned as an output in iopt(IO + 1)%idummy.
Nearest power of 2 ≥ n is returned as an output in iopt(IO + 1)%idummy.
iopt(IO) = ?_options(?_fast_dft_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_dft_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_DFT, D_FAST_DFT, C_FAST_DFT, and Z_FAST_DFT.
Description
The fast_dft routine is a Fortran 90 version of the FFT suite of IMSL (1994, pp. 772-776). The maximum computing efficiency occurs when the size of the array can be factored in the form
using non-negative integer values {i1, i2, i3, i4}. There is no further restriction on n ≥ 1.
Fatal and Terminal Messages
See the messages.gls file for error messages for FAST_DFT. These error messages are numbered 651–661; 701–711.
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_dft_int
use rand_gen_int
implicit none
! This is Example 1 for FAST_DFT.
integer, parameter :: n=1024
real(kind(1e0)), parameter :: one=1e0
real(kind(1e0)) err, y(2*n)
complex(kind(1e0)), dimension(n) :: a, b, c
! Generate a random complex sequence.
call rand_gen(y)
a = cmplx(y(1:n),y(n+1:2*n),kind(one))
c = a
! Transform and then invert the sequence back.
call c_fast_dft(forward_in=a, &
forward_out=b)
call c_fast_dft(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_DFT is correct.'
end if
end
Output
Example 1 for FAST_DFT is correct.
Example 2: Cyclical Data with a Linear Trend
This set of data is sampled from a function x(t) = at + b + y(t), where y(t) is a harmonic series. The independent variable is normalized as –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 LIN_SOL_LSQ, Chapter 1. Then, the residuals are transformed and the resulting frequencies are analyzed.
use fast_dft_int
use lin_sol_lsq_int
use rand_gen_int
use sort_real_int
implicit none
! This is Example 2 for FAST_DFT.
integer i
integer, parameter :: n=64, k=4
integer ip(n)
real(kind(1e0)), parameter :: one=1e0, two=2e0, zero=0e0
real(kind(1e0)) delta_t, pi
real(kind(1e0)) y(k), z(2), indx(k), t(n), temp(n)
complex(kind(1e0)) a_trend(n,2), a, b_trend(n,1), b, c(k), f(n),&
r(n), x(n), x_trend(2,1)
! Generate random data for linear trend and harmonic series.
call rand_gen(z)
a = z(1); b = z(2)
call rand_gen(y)
! This emphasizes harmonics 2 through k+1.
c = y + one
! Determine sampling interval.
delta_t = two/n
t=(/(-one+i*delta_t, i=0,n-1)/)
! Compute pi.
pi = atan(one)*4E0
indx=(/(i*pi,i=1,k)/)
! Make up data set as a linear trend plus harmonics.
x = a + b*t + &
matmul(exp(cmplx(zero,spread(t,2,k)*spread(indx,1,n),kind(one))),c)
! Define least-squares matrix data for a linear trend.
a_trend(1:,1) = one
a_trend(1:,2) = t
b_trend(1:,1) = x
! 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/))
! Transform harmonic residuals.
call c_fast_dft(forward_in=r, forward_out=f)
ip=(/(i,i=1,n)/)
! The dominant frequencies should be 2 through k+1.
! Sort the magnitude of the transform first.
call s_sort_real(-(abs(f)), temp, iperm=ip)
! The dominant frequencies are output in ip(1:k).
! Sort these values to compare with 2 through k+1.
call s_sort_real(real(ip(1:k)), temp)
ip(1:k)=(/(i,i=2,k+1)/)
! Check the results.
if (count(int(temp(1:k)) /= ip(1:k)) == 0) then
write (*,*) 'Example 2 for FAST_DFT is correct.'
end if
end
Output
Example 2 for FAST_DFT is correct.
Example 3: Several 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_dft.
use fast_dft_int
use rand_gen_int
implicit none
! This is Example 3 for FAST_DFT.
! The value of the array size for work(:) is computed in the
! routine fast_dft as a first step.
integer, parameter :: n=64
integer ido_value
real(kind(1e0)) :: one=1e0
real(kind(1e0)) err, y(2*n)
complex(kind(1e0)), dimension(n) :: a, b, save_a
complex(kind(1e0)), allocatable :: work(:)
! Generate a random complex array.
call rand_gen(y)
a = cmplx(y(1:n),y(n+1:2*n),kind(one))
save_a = 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))
call c_fast_dft(forward_in=a, forward_out=b, &
ido=ido_value, work_array=work)
if (ido_value == 1) exit
end do
! Re-enter routine with working values available in work(:).
call c_fast_dft(inverse_in=b, inverse_out=a, &
ido=ido_value, work_array=work)
! Deallocate the space used for work(:).
if (allocated(work)) deallocate(work)
! Check the results.
err = maxval(abs(save_a-a))/maxval(abs(save_a))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 3 for FAST_DFT is correct.'
end if
end
Output
Example 3 for FAST_DFT is correct.
Example 4: Convolutions using Fourier Transforms
In this example we compute sums
The definition implies a matrix-vector product. A direct approach requires about operations consisisting of an add and multiply. An efficient method consisting of computing the products of the transforms of the
{aj} and {bj}
then inverting this product, is preferable to the matrix-vector approach for large problems. The example is also illustrated in operator_ex37, Chapter 10 using the generic function interface FFT and IFFT.
use fast_dft_int
use rand_gen_int
implicit none
! This is Example 4 for FAST_DFT.
integer j
integer, parameter :: n=40
real(kind(1e0)) :: one=1e0
real(kind(1e0)) err
real(kind(1e0)), dimension(n) :: x, y, yy(n,n)
complex(kind(1e0)), dimension(n) :: a, b, c, d, e, f
! Generate two random complex sequence 'a' and 'b'.
call rand_gen(x)
call rand_gen(y)
a=x; b=y
! Compute the convolution 'c' of 'a' and 'b'.
! Use matrix times vector for test results.
yy(1:,1)=y
do j=2,n
yy(2:,j)=yy(1:n-1,j-1)
yy(1,j)=yy(n,j-1)
end do
c=matmul(yy,x)
! Transform the 'a' and 'b' sequences into 'd' and 'e'.
call c_fast_dft(forward_in=a, &
forward_out=d)
call c_fast_dft(forward_in=b, &
forward_out=e)
! Invert the product d*e.
call c_fast_dft(inverse_in=d*e, &
inverse_out=f)
! Check the Convolution Theorem:
! inverse(transform(a)*transform(b)) = convolution(a,b).
err = maxval(abs(c-f))/maxval(abs(c))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 4 for FAST_DFT is correct.'
end if
end
Output
Example 4 for FAST_DFT is correct.
