Fourier Transform (1-D Serial)

!

!  Copyleft: Rupak Mukherjee, Institute for Plasma Research, India, April 2, 2017

!  This is a slow serial 1-dimensional fourier transform program.

!  To run this program in terminal: gfortran fourier_transform.f95; ./a.out

!  To view result: gnuplot> p "Input_Function.dat" u 1:2 w lp, "Input_Function_Recovered.dat" u 1:2 w lp

!  To view Fourier Co-efficients: gnuplot> p "Fourier_Coeff.dat" u 1:2 w lp

!

!==============================================================================================

FUNCTION f(x)

IMPLICIT NONE

real*8, parameter :: pi = 3.14159265358979323846d0

REAL*8 f,x,lamda

!

! YOUR INPUT FUNCTION HERE

!

lamda = 2.0d0

!f= dsin(3.0d0*x)

!f = 7.0d0 + 5.0d0*dsin(2.0d0*x)+ 3.0d0*dcos(8.0d0*x/lamda)

if ((x .ge. 0.0d0) .and. (x .le. 5.0d0)) then

  f = 4.0d0

  else

  f = -1.0d0

endif

END FUNCTION

!==============================================================================================

PROGRAM Fourier_Transform

implicit none

real*8, parameter :: pi = 3.14159265358979323846d0

integer, parameter :: N = 1000 ! Number of data points.

integer, parameter :: UN = N/2+1 ! You may use lower mode-numbers to check convergence.

real*8, parameter :: Length = 6.0d0*pi ! Maximum value of Independent variable.

real*8 dx,f,x0,a0,a,b,k,ft,Re,Im

real*8 x(N),y(0:N),y_rec(0:N),rfft(N),ifft(N)

integer i,u

external f

!open(unit=1,file='Data_File.dat',status='old')

open(unit=10,file='Input_Function.dat',status='unknown')

open(unit=20,file='Fourier_Coeff.dat',status='unknown')

open(unit=30,file='Input_Function_Recovered.dat',status='unknown')

x0 = 0.0d0

dx = Length/dfloat(N)

y(0) = f(x0)

do i = 1,N-1

  x(i) = x0 + dx*dfloat(i)

  y(i) = f(x(i))

!  read(1,*) x(i), y(i)

  write(10,*) x(i), y(i)

enddo ! i

a0 = y(0)

do i = 1,N-1

  a0 = a0 + y(i)

enddo ! i

write(20,*) 0, a0/dfloat(N)

do u = 1,UN-1

  k = 2.0d0*pi*dfloat(u)/Length

  rfft(u) = f(x0)

  ifft(u) = 0.0d0

    do i = 1,N-1

      a = 2.0d0*pi*x(i)/Length

      rfft(u) = rfft(u) + y(i)*dcos(dfloat(u)*a)

      ifft(u) = ifft(u) + y(i)*dsin(dfloat(u)*a)

    enddo ! i

  b = dsqrt(rfft(u)**2 + ifft(u)**2)/dfloat(N)

  ft = 2.0d0*b

  write(20,*) k, ft!, rfft(u)/dfloat(N), ifft(u)/dfloat(N)

enddo ! u

do u = N-1,(N-UN)+1,-1

!  k = 2.0d0*pi*dfloat(u)/Length

  rfft(u) = f(x0)

  ifft(u) = 0.0d0

    do i = 1,N-1

      a = 2.0d0*pi*x(i)/Length

      rfft(u) = rfft(u) + y(i)*dcos(dfloat(u)*a)

      ifft(u) = ifft(u) + y(i)*dsin(dfloat(u)*a)

    enddo ! i

!  b = dsqrt(rfft(u)**2 + ifft(u)**2)/dfloat(N)

!  ft = 2.0d0*b

!  write(20,*) k, ft!, rfft(u)/dfloat(N), ifft(u)/dfloat(N)

enddo ! u

do i = 1,N-1

  Re = 0.0d0

  Im = 0.0d0

    do u = 1,UN-1

      Re = Re + rfft(u)*dcos(2.0d0*pi*dfloat(u)*x(i)/Length)

      Im = Im + ifft(u)*dsin(2.0d0*pi*dfloat(u)*x(i)/Length)

    enddo ! u

    do u = N-1,(N-UN)+1,-1

      Re = Re + rfft(u)*dcos(2.0d0*pi*dfloat(u)*x(i)/Length)

      Im = Im + ifft(u)*dsin(2.0d0*pi*dfloat(u)*x(i)/Length)

    enddo

  y_rec(i) = (a0 + Re + Im)/dfloat(N)

  write(30,*) x(i), y_rec(i)

enddo ! i

END PROGRAM Fourier_Transform