Welcome to AE Resources

Converted document FAST FOURIER TRANSFORMS DEMONSTRATION

FAST FOURIER TRANSFORMS DEMONSTRATION

This code demonstrates the FFT of a set of discrete data. The data set had 361 data points. Since the FFT requires data elements to the power of two, zeros padded the data array from 362 to 512, the next number using the power of 2.

PROGRAM FFT_discrete
!Description:  This program determines the FFT of data
!     provided from a .dat file.

IMPLICIT NONE
REAL, PARAMETER :: pi=4.0*ATAN(1.0)
INTEGER :: i, nl, ierr
REAL :: fs
CHARACTER(LEN=50) :: myfile
INTEGER, DIMENSION(:), ALLOCATABLE :: f
REAL, DIMENSION(:), ALLOCATABLE :: time, myreal, myimag
!COMPLEX, DIMENSION(:,:), ALLOCATABLE :: comp
REAL, DIMENSION(:,:), ALLOCATABLE :: comp
myfile = ’compData.dat’ !Computational data file

CALL readFile(myfile)
WRITE(*,*)

!Allocate
ALLOCATE(time(nl),myreal(nl),myimag(nl), f(nl/2))

!Notes:
!  w = 2*pi*freq
!  f = 1/T   where T = period

!Set time precision
!Implied DO Loop used here!
time(1:nl) = (/((i)*(1/REAL(nl)), i=0,nl-1)/)  !Time for both cos(2*pi*t)
!Set Sampling Frequency
fs = SIZE(time)/time(SIZE(time))

!Split up real and imaginary portions of function
myreal(1:nl) = comp(1:nl,2)
myimag(1:nl) = (/((0.*i), i=1,nl)/)

!Run it through the FFT via scramble and unscramble
CALL scram(512,myreal,myimag)
CALL unscram(512,myreal, myimag)
!Note:: "512" used because it is the next power 
! of 2 higher than the length (361)

!Store magnitudes by reusing myreal vector
myreal(1:nl) = (/(SQRT(myreal(i)**2 + myimag(i)**2), i=1,nl)/)

!Create frequency axis values
f(1:nl/2) = (/(REAL(i)/REAL(SIZE(time)),&
 i=0,SIZE(time)*(int(fs/2.)-1),int(fs))/)
WRITE(*,*) ’FFT Transformation:’
WRITE(*,*) ’(Frequency & Magnitude)’
DO i = 1,nl/2-1
	WRITE(*,*) f(i), myreal(i)
ENDDO


!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! SUBROUTINES !!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!
CONTAINS
SUBROUTINE readFile(name)
	IMPLICIT NONE
	INTEGER :: i, ios
	CHARACTER(LEN=50), INTENT(IN) :: name

	WRITE(*,*) !Extra line for aesthetics
!Read from computed data
	WRITE(*,*) ’Reading from file: ’, name
	!Open COMPUTATIONAL DATA file to read from it
	OPEN(UNIT=2, FILE=name)

	!Determine the number of lines in the file
	nl = 0 !Initialize
	DO i = 1,1000,1
		READ (2, *, IOSTAT=ios) !Dummy variable
		!If it’s the end of the file, then exit
		IF (ios /= 0) THEN
			EXIT
		ELSEIF (i == 1000) THEN
			WRITE (*,*) ’Maximum allowable data (1000 lines)’
		ELSE
			nl = nl + 1
		ENDIF
	ENDDO
	WRITE(*,*) ’   End of file found:’, nl, ’pieces of data total.’
	CLOSE(2)