Module_FFT_Procedures¶
Copyright 2024 IRD
This file is part of statpack.
statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.
Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)
MODULE EXPORTING FAST FOURIER TRANSFORMS.
LATEST REVISION : 26/02/2024
subroutine init_fft ( shap, dim )¶
Purpose¶
Subroutine INIT_FFT sets up constants, the Chirp functions and the Fourier transform of the Chirp functions for use by generic subroutines FFT, FFTXY, FFT_ROW and REAL_FFT for a complex valued array of shape SHAP.
Arguments¶
- SHAP (INPUT) integer(i4b), dimension(:)
- Rank-one integer holding the shape of the complex valued array to be transformed. Size( SHAP ) must be less or equal to 3.
- DIM (INPUT, OPTIONAL) integer(i4b)
- Eventually specifies the index for the Fourier transform. Fourier transform on DIM-index-sections, only. DIM must be less or equal to size( SHAP ).
Further Details¶
INIT_FFT is first called to establish and transform the Chirp functions and other constants. Then, subroutines FFT, FFTXY, FFT_ROW and REAL_FFT can be called any number of times without the precalculated constants being destroyed; a further call to INIT_FFT will only be necessary if Fourier transforms for a new length (or shape) are required.
subroutine init_fft ( length1 )¶
Purpose¶
Subroutine INIT_FFT sets up constants, the bit reverse tables, the Chirp function and the Fourier transform of the Chirp function for use by generic subroutines FFT and FFT_ROW for a series of length LENGTH1.
Arguments¶
- LENGTH1 (INPUT) integer(i4b)
- The length of the complex valued sequence to be transformed. LENGTH1 may be any positive integer.
Further Details¶
INIT_FFT is first called to establish and transform the Chirp function and other constants. Then, subroutine FFT (or FFT_ROW) can be called any number of times without the precalculated constants being destroyed; a further call to INIT_FFT will only be necessary if a new length is required.
This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
subroutine init_fft ( length1, length2 )¶
Purpose¶
Subroutine INIT_FFT sets up constants, the Chirp functions and the Fourier transforms of the Chirp functions for use by generic subroutines FFT or FFTXY for a complex matrix of shape (LENGTH1,LENGTH2).
Arguments¶
- LENGTH1 (INPUT) integer(i4b)
- The number of rows of the complex matrix to be transformed. LENGTH1 may be any positive integer.
- LENGTH2 (INPUT) integer(i4b)
- The number of columns of the complex matrix to be transformed. LENGTH2 may be any positive integer.
Further Details¶
INIT_FFT is first called to establish and transform the Chirp functions and other constants. Then, subroutine FFT (or FFTXY) can be called any number of times without the precalculated constants being destroyed; a further call to INIT_FFT will only be necessary if a new shape is required.
This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
subroutine init_fft ( length1, length2, length3 )¶
Purpose¶
Subroutine INIT_FFT sets up constants, the Chirp functions and the Fourier transforms of the Chirp functions for use by generic subroutines FFT or FFTXY for a complex 3d array of shape (LENGTH1,LENGTH2,LENGTH3).
Arguments¶
- LENGTH1 (INPUT) integer(i4b)
- The extent in the first dimension of the complex array to be transformed. LENGTH1 may be any positive integer.
- LENGTH2 (INPUT) integer(i4b)
- The extent in the second dimension of the complex array to be transformed. LENGTH2 may be any positive integer.
- LENGTH3 (INPUT) integer(i4b)
- The extent in the third dimension of the complex array to be transformed. LENGTH3 may be any positive integer.
Further Details¶
INIT_FFT is first called to establish and transform the Chirp functions and other constants. Then, subroutine FFT (or FFTXY) can be called any number of times without the precalculated constants being destroyed; a further call to INIT_FFT will only be necessary if a new shape is required.
This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
subroutine real_fft ( vec, vect, forward)¶
Purpose¶
Subroutine REAL_FFT computes the Fast Fourier Transform (FFT) for a real valued sequence VEC of even length.
Arguments¶
- VEC (INPUT) real(stnd), dimension(:)
On entry, the real valued sequence to be transformed.
Size( VEC ) must be an even (positive) integer.
- VECT (OUTPUT) complex(stnd), dimension(:)
On exit, a complex vector of length size(VEC)/2+1 containing the first size(VEC)/2+1 coefficients of the Fourier transform of the real sequence VEC. These coefficients are the positive frequency half of the full complex Fourier transform of the the real value sequence VEC.
VECT must verify: size( VECT ) = size( VEC )/2 + 1.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
REAL_FFT computes the Forward Fourier Transform, VECT, according to the following formula
VECT(j+1) = [ sum k=0 to nn-1 ] VEC(k+1) exp( - i 2 pi j k / nn )for j=0, 1, …, nn/2 and where i=sqrt( -1 ), nn=size(VEC) and pi=3.1415923565…
REAL_FFT computes the Backward Fourier Transform, VECT, according to the following formula
VECT(j+1) = (1/nn) [ sum k=0 to nn-1 ] VEC(k+1) exp( 2 pi j k / nn )for j=0, 1, …, nn/2 and where i=sqrt( -1 ), nn=size(VEC) and pi=3.1415923565…
The remaining values of the Fourier Transform may be computed by the following lines of code
nn = size(VEC)
nnd2 = nn/2
vect(nn:nnd2+2:-1) = conjg( vect(2:nnd2) )
Before using REAL_FFT, the user must call subroutine INIT_FFT as follows :
call init_fft( size(vec)/2 )For more details on the Discrete Fourier Transform, see:
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
- Bloomfield, P., 1976:
- Fourier analysis of time series- An introduction. John Wiley and Sons, New York.
subroutine real_fft ( mat, matt, forward)¶
Purpose¶
Subroutine REAL_FFT computes the Fast Fourier Transform (FFT) for each row of a real valued matrix MAT.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, the real valued matrix to be transformed. Size( MAT , 2 ) must be an even (positive) integer.
- MATT (OUTPUT) complex(stnd), dimension(:,:)
On exit, a complex matrix of shape size(MAT,1) by size(MAT,2)/2+1. each row of MATT contains the first size(MAT,2)/2+1 coefficients of the Fourier transform of the corresponding row of the real matrix MAT. These coefficients are the positive frequency half of the full complex Fourier transform of the corresponding row of MAT.
The shape of MATT must verify:
- size( MATT, 1 ) = size( MAT, 1 )
- size( MATT, 2 ) = size( MAT, 2 )/2 + 1.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
REAL_FFT computes the Forward Fourier Transform, MATT, according to the following formula
MATT(l,j+1) = [ sum k=0 to nn-1 ] MAT(l,k+1) exp( - i 2 pi j k / nn )for j=0, 1, …, nn/2 , l=1, …, size(MAT,1) and where i=sqrt( -1 ), nn=size(VEC) and pi=3.1415923565…
REAL_FFT computes the Backward Fourier Transform, MATT, according to the following formula
MATT(l,j+1) = (1/nn) [ sum k=0 to nn-1 ] MAT(l,k+1) exp( 2 pi j k / nn )for j=0, 1, …, nn/2 , l=1, …, size(MAT,1) and where i=sqrt( -1 ), nn=size(VEC) and pi=3.1415923565…
The remaining values of the Fourier Transform may be computed by the following lines of code
nn = size(mat,2)
nnd2 = nn/2
matt(l,nn:nnd2+2:-1) = conjg( matt(l,2:nnd2) )
Before using REAL_FFT, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(MAT,1), size(MAT,2)/2 /), dim=2)For more details on the Discrete Fourier Transform, see:
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
- Bloomfield, P., 1976:
- Fourier analysis of time series- An introduction. John Wiley and Sons, New York.
subroutine real_fft_forward ( vec, vecr, veci)¶
Purpose¶
Subroutine REAL_FFT_FORWARD implements the forward discrete Fourier Transform for a real valued sequence VEC of general length.
Arguments¶
- VEC (INPUT) real(stnd), dimension(:)
- On entry, the real valued sequence to be transformed.
- VECR (OUTPUT) real(stnd), dimension(:)
On exit, the real part of the forward discrete Fourier Transform of the sequence VEC.
VECR must verify: size( VECR ) = size( VEC )/2 + 1.
- VECI (OUTPUT) real(stnd), dimension(:)
On exit, the imaginary part of the forward discrete Fourier Transform of the sequence VEC.
VECI must verify: size( VECI ) = size( VEC )/2 + 1.
Further Details¶
Only, the part of the discrete Fourier Transform corresponding to the positive frequencies are computed and output in VECR and VECI.
The forward Discrete Fourier Transform is computed using Goertzel method.
For more details on the Goertzel method for computing the Discrete Fourier Transform. For further details, see:
- Goertzel, G., 1958:
- An Algorithm for the Evaluation of Finite Trigonometric Series, The American Mathematical Monthly, Vol. 65, No. 1, pp. 34-35
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine real_fft_backward ( vecr, veci, vec)¶
Purpose¶
Subroutine REAL_FFT_BACKWARD implements the (real) backward discrete Fourier Transform for a complex valued sequence stored in the vector VECR (real part of the complex sequence) and VECI (imaginary part of the sequence). The resulting real discrete Fourier Transform is stored in the real vector VEC.
Size(VEC) which gives the size of the transform may be of general length.
Arguments¶
- VECR (INPUT) real(stnd), dimension(:)
On entry, the real part of the complex sequence to be transformed.
VECR must verify: size( VECR ) = size( VEC )/2 + 1.
- VECI (INPUT) real(stnd), dimension(:)
On entry, the imaginary part of the complex sequence to be transformed.
VECI must verify: size( VECI ) = size( VEC )/2 + 1.
- VEC (OUTPUT) real(stnd), dimension(:)
- On exit, the discrete Fourier transform real valued sequence.
Further Details¶
The backward Discrete Fourier Transform is computed using Goertzel method.
For more details on the Goertzel method for computing the Discrete Fourier Transform. For further details, see:
- Goertzel, G., 1958:
- An Algorithm for the Evaluation of Finite Trigonometric Series, The American Mathematical Monthly, Vol. 65, No. 1, pp. 34-35
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine real_fft_forward ( mat, matr, mati, dim )¶
Purpose¶
Subroutine REAL_FFT_FORWARD computes the forward discrete Fourier Transform of each row (DIM=2) or each column (DIM=1) of the real matrix MAT.
Size(MAT,DIM) which gives the size of the transform may be of general length.
The real parts of the forward discrete Fourier Transforms are stored (rowwise) in MATR and the corresponding imaginary parts of the transforms are stored in MATI.
Arguments¶
- MAT (INPUT) real(stnd), dimension(:,:)
- On entry, the real valued sequences to be transformed.
- MATR (OUTPUT) real(stnd), dimension(:,:)
On exit, the real part of the forward discrete Fourier Transform of the sequences stored in MAT.
The shape of MATR must verify:
- size( MATR, 1 ) = size( MAT, 3-DIM )
- size( MATR, 2 ) = size( MAT, DIM )/2 + 1.
- MATI (OUTPUT) real(stnd), dimension(:,:)
On exit, the imaginary part of the forward discrete Fourier Transform of the sequences stored in MAT.
The shape of MATI must verify:
- size( MATI, 1 ) = size( MAT, 3-DIM )
- size( MATI, 2 ) = size( MAT, DIM )/2 + 1.
- DIM (INPUT) integer(i4b)
Specifies the index for the Fourier transform. If:
- DIM = 1 : Fourier transform on first index.
- DIM = 2 : Fourier transform on second index.
Further Details¶
Only, the parts of the discrete Fourier Transforms corresponding to the positive frequencies are computed and output in MATR and MATI.
The forward Discrete Fourier Transform is computed using Goertzel method.
For more details on the Goertzel method for computing the Discrete Fourier Transform. For further details, see:
- Goertzel, G., 1958:
- An Algorithm for the Evaluation of Finite Trigonometric Series, The American Mathematical Monthly, Vol. 65, No. 1, pp. 34-35
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine real_fft_backward ( matr, mati, mat, dim )¶
Purpose¶
Subroutine REAL_FFT_BACKWARD computes the (real) backward discrete Fourier Transform for complex valued sequences stored in the matrices MATR (real part of the sequences stored rowwise) and MATI (imaginary part of the sequences stored rowwise). The resulting real discrete Fourier Transforms are stored in the rows (DIM=2) or the columns (DIM=1) of the real matrix MAT.
Size(MAT,DIM) which gives the size of the transform may be of general length.
Arguments¶
- MATR (INPUT) real(stnd), dimension(:,:)
On entry, the real part of the complex sequences to be transformed.
The shape of MATR must verify:
- size( MATR, 1 ) = size( MAT, 3-DIM )
- size( MATR, 2 ) = size( MAT, DIM )/2 + 1.
- MATI (INPUT) real(stnd), dimension(:,:)
On entry, the imaginary part of the complex sequences to be transformed.
The shape of MATI must verify:
- size( MATI, 1 ) = size( MAT, 3-DIM )
- size( MATI, 2 ) = size( MAT, DIM )/2 + 1.
- MAT (OUTPUT) real(stnd), dimension(:,:)
- On exit, the real backward discrete Fourier transforms of the complex sequences stored rowwise in MATR and MATI.
- DIM (INPUT) integer(i4b)
Specifies the index for the Fourier transform. If:
- DIM = 1 : Fourier transform on first index.
- DIM = 2 : Fourier transform on second index.
Further Details¶
The backward Discrete Fourier Transform is computed using Goertzel method.
For more details on the Goertzel method for computing the Discrete Fourier Transform. For further details, see:
- Goertzel, G., 1958:
- An Algorithm for the Evaluation of Finite Trigonometric Series, The American Mathematical Monthly, Vol. 65, No. 1, pp. 34-35
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine fftxy ( x, y, fftx, ffty)¶
Purpose¶
Given two real valued sequences of the same length, X and Y, FFTXY returns the Fast Fourier Transforms of these sequences in the two complex valued sequences FFTX and FFTY.
Arguments¶
- X, Y (INPUT) real(stnd), dimension(:)
On entry, the real valued sequences to be transformed.
X and Y must verify: size( X ) = size( Y ).
- FFTX, FFTY (OUTPUT) complex(stnd), dimension(:)
On exit, FFTX, FFTY are replaced by the Fourier transforms of X and Y, respectively.
FFTX and FFTY must verify: size( FFTX ) = size( FFTY ) = size( X ) = size( Y ).
Further Details¶
Size( FFTX ) = size( FFTY ) = size( X ) = size( Y ) may be of general length.
If size(X) is an exact power of two, Bailey’s Four-Step FFT algorithm is used, otherwise the CHIRP-Z transform is employed.
Before using FFTXY, the user must call subroutine INIT_FFT as follows :
call init_fft( size(X) )This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83 and Bailey (1990). For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Bailey, D., 1990:
- FFTs in External or Hierarchical Memory, The Journal of Supercomputing, 4, 23-35.
subroutine fftxy ( x, y, fftx, ffty )¶
Purpose¶
Given two real valued matrices of the same shape, X and Y, FFTXY returns the Fast Fourier Transforms of X and Y in the two complex valued matrices FFTX and FFTY, respectively.
Arguments¶
- X, Y (INPUT) real(stnd), dimension(:,:)
On entry, the real valued matrices to be transformed.
X and Y must verify the equality: shape( X ) = shape( Y ).
- FFTX, FFTY (OUTPUT) complex(stnd), dimension(:,:)
On exit, FFTX and FFTY are replaced by the Fourier transforms of X and Y, respectively.
FFTX and FFTY must verify the equalities:
- shape( FFTX ) = shape( FFTY ) = shape( X ) = shape( Y ).
Further Details¶
Depending if size(X,1) and size(X,2) are exact powers of two or not, a radix-2 decimation-in-time Cooley-Tukey algorithm or a CHIRP-Z transform is employed.
Before using FFTXY, the user must call subroutine INIT_FFT as follows :
call init_fft( size(X,1), size(X,2) )For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Cooley, J.W., Lewis, P., and Welch, P., 1969:
- The Fast Fourier Transform and its Applications. IEEE Trans on Education, 12, 1, 28-34.
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine fftxy ( x, y, fftx, ffty )¶
Purpose¶
Given two real valued 3D arrays of the same shape, X and Y, FFTXY returns the Fast Fourier Transforms of X and Y in the two complex valued 3D arrays FFTX and FFTY, respectively.
Arguments¶
- X, Y (INPUT) real(stnd), dimension(:,:,:)
On entry, the real valued 3D arrays to be transformed.
X and Y must verify: shape( X ) = shape( Y ).
- FFTX, FFTY (OUTPUT) complex(stnd), dimension(:,:,:)
On exit, FFTX and FFTY are replaced by the Fourier transforms of X and Y, respectively.
FFTX and FFTY must verify: shape( FFTX ) = shape( FFTY ) = shape( X ) = shape( Y ).
Further Details¶
Depending if size(X,1), size(X,2) and size(X,3) are exact powers of two or not, a radix-2 decimation-in-time Cooley-Tukey algorithm or a CHIRP-Z transform is employed.
Before using FFTXY, the user must call subroutine INIT_FFT as follows :
call init_fft( size(X,1), size(X,2), size(X,3) )For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Cooley, J.W., Lewis, P., and Welch, P., 1969:
- The Fast Fourier Transform and its Applications. IEEE Trans on Education, 12, 1, 28-34.
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine fftxy ( x, y, fftx, ffty, dim )¶
Purpose¶
Given two real valued matrices of the same shape, X and Y, FFTXY returns the Fast Fourier Transforms of the rows (DIM=2) or the columns (DIM=1) of X and Y in the two complex valued matrices FFTX and FFTY, respectively.
Arguments¶
- X, Y (INPUT) real(stnd), dimension(:,:)
On entry, the real valued matrices to be transformed.
X and Y must verify: shape( X ) = shape( Y ).
- FFTX, FFTY (OUTPUT) complex(stnd), dimension(:,:)
On exit:
- each row of FFTX and FFTY are replaced by the Fourier transforms of the rows of X and Y, respectively, if DIM=2 .
- each column of FFTX and FFTY are replaced by the Fourier transforms of the colums of X and Y, respectively, if DIM=1 .
FFTX and FFTY must verify: shape( FFTX ) = shape( FFTY ) = shape( X ) = shape( Y ).
- DIM (INPUT) integer(i4b)
Specifies the index for the Fourier transform. If:
- DIM = 1 : Fourier transform on first index,
- DIM = 2 : Fourier transform on second index.
Further Details¶
Size(FFTX,DIM) = size(FFTY,DIM) = size(X,DIM) = size(Y,DIM) may be of general length.
Before using FFTXY, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(X,1), size(X,2) /), dim=DIM )
subroutine fftxy ( x, y, fftx, ffty, dim )¶
Purpose¶
Given two real valued 3D arrays of the same shape, X and Y, FFTXY returns the Fast Fourier Transforms of each DIM-index section of X and Y in the two complex valued 3D arrays FFTX and FFTY, respectively.
Arguments¶
- X, Y (INPUT) real(stnd), dimension(:,:,:)
On entry, the real valued 3D arrays to be transformed.
X and Y must verify: shape( X ) = shape( Y ).
- FFTX, FFTY (OUTPUT) complex(stnd), dimension(:,:,:)
On exit, the DIM-index sections of FFTX and FFTY are replaced by the Fourier transforms of the DIM-index sections of X and Y, respectively.
FFTX and FFTY must verify: shape( FFTX ) = shape( FFTY ) = shape( X ) = shape( Y ).
- DIM (INPUT) integer(i4b)
Specifies the index for the Fourier transform. If:
- DIM = 1 : Fourier transform on first-index-sections.
- DIM = 2 : Fourier transform on second-index-sections.
- DIM = 3 : Fourier transform on third-index-sections.
Further Details¶
Size(FFTX,DIM) = size(FFTY,DIM) = size(X,DIM) = size(Y,DIM) may be of general length.
Before using FFTXY, the user must call subroutine INIT_FFT as follows:
call init_fft( (/ size(X,1), size(X,2), size(X,3) /), dim=DIM )
subroutine fft ( dat, forward)¶
Purpose¶
Subroutine FFT implements the Fast Fourier Transform for a complex valued sequence DAT of general length.
Forward discrete Fourier transform of a vector DAT(:) is given by
t( DAT )(j) = [ sum k=0 to nn-1 ] DAT(k) exp( - i 2 pi j k / nn )
Backward discrete Fourier transform of a vector DAT(:) is given by
t( DAT )(j) = (1/nn) [ sum k=0 to nn-1 ] DAT(k) exp( i 2 pi j k / nn )
where i = sqrt( -1 ), nn = size(DAT) and pi = 3.1415923565…
Note that the indexing of DAT is shifted by one : DAT(0) stored in DAT(1), … , DAT( nn-1 ) stored in DAT( nn ).
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:)
- On entry, the complex valued sequence to be transformed. On exit, DAT is replaced by the Fourier transform.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
If size(DAT) is an exact power of two, Bailey’s Four-Step FFT algorithm is used, otherwise the CHIRP-Z transform is employed.
Before using FFT, the user must call subroutine INIT_FFT as follows :
call init_fft( size(DAT) )This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83 and Bailey (1990). For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Bailey, D., 1990:
- FFTs in External or Hierarchical Memory. The Journal of Supercomputing, 4, 23-35.
subroutine fft ( dat, forward)¶
Purpose¶
Subroutine FFT implements the Fast Fourier Transform for a complex matrix DAT of general shape.
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:,:)
- On entry, the complex matrix to be transformed. On exit, DAT is replaced by its Fourier transform.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
Depending if size(DAT,1) and size(DAT,2) are exact powers of two or not, a radix-2 decimation-in-time Cooley-Tukey algorithm or a CHIRP-Z transform is employed.
Before using FFT, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(DAT,1), size(DAT,2) /) )For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Cooley, J.W., Lewis, P., and Welch, P., 1969:
- The Fast Fourier Transform and its Applications. IEEE Trans on Education, 12, 1, 28-34.
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine fft ( dat, forward)¶
Purpose¶
Subroutine FFT implements the Fast Fourier Transform for a complex 3D array DAT of general shape.
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:,:,:)
- On entry, the complex array to be transformed. On exit, DAT is replaced by its Fourier transform.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
Depending if size(DAT,1), size(DAT,2) and size(DAT,3) are exact powers of two or not, a radix-2 decimation-in-time Cooley-Tukey algorithm or a CHIRP-Z transform is employed.
Before using FFT, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(DAT,1), size(DAT,2), size(DAT,3) /) )For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Cooley, J.W., Lewis, P., and Welch, P., 1969:
- The Fast Fourier Transform and its Applications. IEEE Trans on Education, 12, 1, 28-34.
- Oppenheim, A.V., and Schafer, R.W., 1999:
- Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.
subroutine fft ( dat, forward, dim)¶
Purpose¶
Subroutine FFT replaces each row of DAT by its Fourier transform. (DIM=2) or each column of DAT by its Fourier transform (DIM=1). Size(DAT,DIM) may be of general length.
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:,:)
- On entry, the complex valued sequences to be transformed. On exit, each row of DAT is replaced by its Fourier transform if DIM=2 or each column of DAT is replaced by its Fourier transform if DIM=1.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
- DIM (INPUT) integer(i4b)
Specifies the index for the Fourier transform. If:
- DIM = 1 : Fourier transform on first index,
- DIM = 2 : Fourier transform on second index.
Further Details¶
If size(DAT,DIM) is an exact power of two, a 1D in-place complex-complex radix-2 decimation-in-time Cooley-Tukey FFT algorithm is used, otherwise the CHIRP-Z transform is employed.
Before using FFT, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(DAT,1), size(DAT,2) /), dim=DIM )This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
subroutine fft ( dat, forward, dim)¶
Purpose¶
Subroutine FFT replaces each DIM-index section of DAT by its Fourier transform. Size(DAT,DIM) may be of general length.
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:,:,:)
- On entry, the complex valued sequences to be transformed. On exit, the DIM-index sections of DAT are replaced by their Fourier transforms.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
- DIM (INPUT) integer(i4b)
Specifies the index for the Fourier transform. If:
- DIM = 1 : Fourier transform on first-index-sections,
- DIM = 2 : Fourier transform on second-index-sections,
- DIM = 3 : Fourier transform on third-index-sections.
Further Details¶
If size(DAT,DIM) is an exact power of two, a 1D in-place complex-complex radix-2 decimation-in-time Cooley-Tukey FFT algorithm is used, otherwise the CHIRP-Z transform is employed.
Before using FFT_DIM_CT, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(DAT,1), size(DAT,2), size(DAT,3) /), dim=DIM )This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
subroutine fft_row ( dat, forward)¶
Purpose¶
Subroutine FFT_ROW implements the Fast Fourier Transform for a complex valued sequence DAT of general length.
Forward discrete Fourier transform of a vector DAT(:) is given by
t( DAT )(j) = [ sum k=0 to nn-1 ] DAT(k) exp( - i 2 pi j k / nn )
Backward discrete Fourier transform of a vector DAT(:) is given by
t( DAT )(j) = (1/nn) [ sum k=0 to nn-1 ] DAT(k) exp( i 2 pi j k / nn )
where i = sqrt( -1 ), nn = size(DAT) and pi = 3.1415923565…
Note that the indexing of DAT is shifted by one : DAT(0) stored in DAT(1), … , DAT( nn-1 ) stored in DAT( nn ).
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:)
- On entry, the complex valued sequence to be transformed. On exit, DAT is replaced by the Fourier transform.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
If size(DAT) is an exact power of two, Bailey’s Four-Step FFT algorithm is used, otherwise the CHIRP-Z transform is employed.
This is the parallelized version of FFT subroutine (Parallelization is done with OPENMP directives).
Before using FFT_ROW, the user must call subroutine INIT_FFT as follows :
call init_fft( size(DAT) )This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83. For more details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.
- Bailey, D., 1990:
- FFTs in External or Hierarchical Memory. The Journal of Supercomputing, 4, 23-35.
subroutine fft_row ( dat, forward)¶
Purpose¶
Subroutine FFT_ROW replaces each row of DAT by its Fourier transform. Size(DAT,2) may be of general length.
Arguments¶
- DAT (INPUT/OUTPUT) complex(stnd), dimension(:,:)
- On entry, the complex valued sequences to be transformed. On exit, each row of DAT is replaced by its Fourier transform.
- FORWARD (INPUT) logical(lgl)
Specifies whether a forward or backward Fourier transform is desired. If:
- FORWARD = true: a forward Fourier transform is computed
- FORWARD = false: a backward Fourier transform is computed.
Further Details¶
If size(DAT,2) is an exact power of two, a 1D complex-complex radix-2 decimation-in-time Cooley-Tukey FFT algorithm is used, otherwise the CHIRP-Z transform is employed.
This is a parallelized FFT subroutine if OPENMP is used (Parallelization is done with OPENMP directives).
Before using FFT_ROW, the user must call subroutine INIT_FFT as follows :
call init_fft( (/ size(DAT,1), size(DAT,2) /), dim=2 )This subroutine is adapted from Applied Statistics algorithms AS 117 and AS 83. For further details, see:
- Monro, D.M., and Branch, J.L., 1977:
- The Chirp discrete Fourier transform of general length. Appl. Statist., 26 (3), 351-361.