C ------------------------------------------------------------------------- C | File name: fftd.for Sean O'Connor 9-22-87 | C | 7-01-17 | C ------------------------------------------------------------------------- C ------------------------------------------------------------------------- C | VAX 11/785 VMS 4.5 VAX EXTENDED FORTRAN 77 | C | MacBook Pro Ubuntu LTS 16.04 Linux gfortran | C ------------------------------------------------------------------------- C **************************************************************************** C * * C * FFTD * C * * C **************************************************************************** C C C DESCRIPTION C C This program tests the subroutine FFT in the file FFT.FOR It does a C forward discrete Fourier transform (DFT). It then does an inverse C transform of the forward transform to check the roundoff error. C C C INPUT C C The number of complex points in the transform. C Whether to do a forward (1) or inverse (-1) transform. C The data points (complex numbers). C C C OUTPUT C C The DFT of the input. C The inverse DFT of the first DFT. C C C EXAMPLE C C $ ./FFTD C C Enter number of complex points in transform: 4 C C Enter 1 (forward) or -1 (inverse) transform: 1 C C Enter point no. 0 (Real, Imag) > (1,1) C Enter point no. 1 (Real, Imag) > (2,2) C Enter point no. 2 (Real, Imag) > (3,3) C Enter point no. 3 (Real, Imag) > (4,4) C C num_points = 4 C direction = 1 C C Point No. 0 = ( 1.0000000000000 1.0000000000000 ) C Point No. 1 = ( 2.0000000000000 2.0000000000000 ) C Point No. 2 = ( 3.0000000000000 3.0000000000000 ) C Point No. 3 = ( 4.0000000000000 4.0000000000000 ) C C The transform of the input is . . . C C Point No. 0= ( 5.0000000000000 5.0000000000000 ) C Point No. 1= ( -2.0000000000000 0.00000000000000E+000) C Point No. 2= ( -1.0000000000000 -1.0000000000000 ) C Point No. 3= (-0.59604644775391E-007 -2.0000000000000 ) C C The inverse transform of the transform above is . . . C C Point No. 0= ( 1.0000000000000 1.0000000000000 ) C Point No. 1= ( 2.0000000000000 2.0000000000000 ) C Point No. 2= ( 3.0000000000000 3.0000000000000 ) C Point No. 3= ( 4.0000000000000 4.0000000000000 ) C C C INSTALLATION C C To compile, type C C $ gfortran FFTD.FOR FORTRAN FFT.FOR -o FFT C C To run, type C C $ ./FFTD C C C **************************************************************************** C * * C **************************************************************************** PROGRAM fftd C ------------------------------------------------------------------------- C | MACRO DEFINITIONS | C ------------------------------------------------------------------------- PARAMETER ( MAX_NUM_POINTS = 255 ) ! Maximum number of points in FFT. C ------------------------------------------------------------------------- C | LOCAL VARIABLES | C ------------------------------------------------------------------------- INTEGER i, ! Loop counter $ num_points, ! Number of points in the transform $ direction ! Direction of transform 1 for forward, ! -1 for inverse COMPLEX X( 0 : MAX_NUM_POINTS ) ! Array of complex data points ! to be transformed. C ------------------------------------------------------------------------- C | PROGRAM CODE | C ------------------------------------------------------------------------- C C Read the number of points in the FFT C WRITE( 6, 100 ) 100 FORMAT ( ' Enter number of complex points in transform: ' ) READ ( 5 , * ) num_points ! Free format input C C Decide whether to do the forward or inverse FFT C WRITE( 6, 200 ) 200 FORMAT ( //, ' Enter 1 for a forward transform ', /, $ ' or -1 for an inverse transform: ' ) READ ( 5, * ) direction ! Free format input C C Check to see that the number of points in the FFT C is valid. If not, halt. C IF ( (num_points .LE. 0) .OR. $ (num_points .GT. MAX_NUM_POINTS) ) THEN WRITE( 6, 300 ) 300 FORMAT( ' The number of points is out of range ') STOP END IF C C Read in the complex data points, one by one C DO i = 0, num_points - 1 WRITE( 6, 400 ) i 400 FORMAT ( ' Enter point no. ', I3,' (Real, Imag) > ') READ ( 5 , * ) X( i ) ! Free format input END DO C C Print out the points just read in. C WRITE( 6, 500 ) num_points 500 FORMAT ( //, ' num_points = ', I10, / ) WRITE( 6, 600 ) direction 600 FORMAT ( ' direction = ', I10, // ) DO i = 0, num_points - 1 WRITE( 6 , 700 ) i , X ( i ) 700 FORMAT ( ' Point no. ', I3, ' = ( ', 2G22.14E3, ') ' ) END DO C C Call the FFT subroutine. C CALL FFT( num_points, X, direction ) C C Print the transform results. C WRITE( 6, 800 ) 800 FORMAT ( //, ' The transform of the input data is . . . ', // ) DO i = 0, num_points - 1 WRITE( 6, 900 ) i , X( i ) 900 FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' ) END DO C C Call the FFT subroutine with -direction to do the inverse transform. C CALL FFT( num_points, X, -direction ) C C Print the inverse transform coefficients. C WRITE ( 6, 1000 ) 1000 FORMAT ( //,' The inverse transform ', $ 'of the transform above is . . .', // ) DO i = 0, num_points - 1 WRITE ( 6, 1100 ) i , X( i ) 1100 FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' ) END DO C C Test the error checking in the FFT routine. C C Number of points is too large. Currently, MAX_TABLE_SIZE in C subroutine FFT is 20, so we can transform up to 2 ** 21 points. PRINT *, ' Calling FFT with too many points.' CALL FFT( 2 ** 22, X, 1 ) C Number of points is less than 2. PRINT *, ' Calling FFT with 1 point.' CALL FFT( 1, X, 1 ) C Calling FFT with 2^ 15 + 1 points (not a power of 2) PRINT *, ' Calling FFT with 2 ^ 15 + 1 points' CALL FFT( 2 * 15 + 1, X, 1 ) C Direction is not 1 or -1. PRINT *, ' Calling FFT with direction = 2' CALL FFT( 16, X, 2) STOP END ! of program fftd