! Fourier analysis routines - compute the FFT of a signal
! this program does forward and backward transforms
! for forward transforms: assume that the data is a real function
! for backward transforms: assume that the final result will be a real function,
!     even though the data to be transformed will have both real and imaginary parts
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program fourier
   library "sgfunc*","sglib*"
   option nolet
   dim t(1024),y(2048)     ! could make this larger to handle larger data sets
   call readin(t,y,n_points,store$,isign)   ! read in data
   call fft(y,n_points,isign)               ! perform fft
   call calc_freq(t,n_points,isign)         ! calculate corresponding frequencies
   call display(t,y,n_points,isign)         ! display results
   call save(t,y,n_points,store$,isign)     ! save them too
end
! calculate the frequencies that go with the transformed signal
! the time array t() matches the array containing the transform
! so t(i) with i odd is for the cosine transform while i even is for
! the sine transform
sub calc_freq(t(),n,isign)
   if isign = 1 then        ! must be careful to handle forward and backward 
      dt = t(3) - t(1)      ! transforms slightly differently
      for i = 1 to n-1 step 2          ! isign = +1 is a forward transform
         t(i) = (i-1) / (2 * n * dt)
         t(i+1) = t(i)
      next i
   else                                ! isign = -1 is a backward transform
      dt = t(3) - t(1) 
      for i = 1 to n
         t(i) = (i-1) / (n * dt)
      next i
   end if
end sub
! save the transform results
! for forward transforms (isign = +1):
!    the cosine transform goes into the file store$.fft.r 
!    the sine transform into store$.fft.i
!    All of the transform data also go into the file store$.fft with real and imaginary
!    parts alternated - this makes it easy to perform the back transform
!    with the same program
! for backward transforms (isign = -1):
!    store the result in the file store$.bfft
!    note that in this case we assume that the signal is a real function
sub save(t(),y(),n,store$,isign)
   if isign = +1 then                 ! forward transform
      open #1: name store$ & ".fft", organization text, create new
      for i = 1 to n-1 step 2
        print #1: t(i),",",y(i)  
        print #1: t(i+1),",",y(i+1)  
      next i
      close #1
      open #1: name store$ & ".fft.r", organization text, create new
      open #2: name store$ & ".fft.i", organization text, create new
      for i = 1 to n-1 step 2
         print #1: t(i),",",y(i)
         print #2: t(i+1),",",y(i+1)
      next i
      close #1
      close #2
   else                                ! backward transform
      open #1: name store$ & ".bfft", organization text, create new
      for i = 1 to n
        print #1: t(i),",",y(i)  
      next i
      close #1
   end if
end sub
! read the data from a file and put it into arrays t() (time data)
! and y() (the signal data)
! note that for a forward transform of a real data set,
! the signal should be in the file as follows
! t(1), y(1)
! t(2), y(2)
! t(3), y(3)
! etc.
!
! for a back transform, the data will in general have both real (cosine)
! and imaginary (sine) parts, so the data must be store differently
! in this case the real and imaginary parts are to be placed on
! alternating lines
! t(1), y_real(1)
! t(1), y_imag(1)
! t(2), y_real(2)
! t(2), y_imag(2)
! etc.
sub readin(t(),y(),n_points,store$,isign)
   input prompt "name of input file -> ": store$
   input prompt "forward (+1) or backward (-1) transform -> ": isign
   open #1: name store$, organization text, create old
   if isign = 1 then
      i = 1
      n_points = 0
      do
         ask #1: pointer place$
         if place$ = "END" then exit do
         input #1: t(i),y(i)
         y(i+1) = 0               ! fill the imaginary parts with zeros
         t(i+1) = t(i)
         i = i + 2
         n_points = n_points + 1
      loop
      close #1
      k = 1   ! pad with zeros up to the next power of 2 for n_points
      do 
         if 2^k >= n_points then exit do
         k = k + 1
      loop
      for j = i to 2^k
         y(j) = 0
      next j
      n_points = 2^k
      mat redim t(2*n_points),y(2*n_points)          ! trim arrays to size used
   else  ! for the inverse transform assume that the number of data
         ! points is already 2^n where n is an integer
         ! so no padding will be required
         ! data for the inverse transform is stored differently - see above
      i = 1
      n_points = 0
      do
         ask #1: pointer place$
         if place$ = "END" then exit do
         input #1: t(i),y(i)              ! real part
         input #1: t(i+1),y(i+1)          ! imaginary part
         i = i + 2
         n_points = n_points + 2
      loop
      close #1
      mat redim t(2*n_points),y(2*n_points)
   end if
end sub
! display the sine and cosine transforms separately if this is a forward
! transform
sub display(t(),d(),n,isign)
   if isign = 1 then
      dim ycos(0),ysin(0),tcos(0)
      mat redim ycos(n/2),ysin(n/2),tcos(n/2)
      for i = 1 to n/2 ! - 1
         ycos(i) = d(2*i-1)
         ysin(i) = d(2*i)
         tcos(i) = t(2*i-1)
      next i
      call datagraph(tcos,ycos,3,1,"black")    ! graph real part
      call adddatagraph(tcos,ysin,4,2,"red")   ! add imaginary part to graph
   else      ! for an inverse transform plot only the real part
      mat redim t(n),d(n)              ! only have to plot the real part
      call datagraph(t,d,3,1,"black")  ! which is already in d()
   end if
   get key z
end sub
! fourier transform of a complex data set
! n_points = number of data points - n_points must be a power of 2
! f_type = +1 -> forward transform
! f_type = -1 -> backward transform
! the data comes in array y()
! real part and imaginary parts are interleaved
! real part in y(0), y 
! imaginary part of data comes in yi() 
! the result is returned in the arrays yr() and yi()
sub fft(y(),n_points,f_type)
   dim yr(0),yi(0)
   option angle radians
   declare def p_val
   j = 1
   n_p = n_points
   mat redim yr(n_points),yi(n_points)
   for i = 1 to n_points                ! split the data up into real and imaginary
      yr(i) = y(j)                      ! parts 
      yi(i) = y(j+1)
      j = j + 2
   next i
   if f_type = -1 then  ! a forward transform
      j = n_points - 1
      for i = n_points + 1 to n_p
         yr(i) = 0  
         yi(i) = 0 
         j = j - 1
      next i
   end if
   n_power = 1
   do    ! first determine n_power where  n_points = 2^n_power   
      if n_p / 2 <= 1 then exit do
      n_p = n_p / 2
      n_power = n_power + 1
   loop
   n1 = n_power - 1
   n2 = n_points / 2
   for i = 1 to n_power
      k = 0
      do
         for j = 1 to n2
            ex = 2 * pi * p_val(k/(2^n1),n_power) / n_points
            cos_factor = cos(ex)
            sin_factor = sin(-f_type*ex)
            k = k + 1
            tmp_real = cos_factor * yr(k+n2) + sin_factor * yi(k+n2)
            tmp_imag = cos_factor * yi(k+n2) - sin_factor * yr(k+n2)
            yr(k+n2) = yr(k) - tmp_real
            yi(k+n2) = yi(k) - tmp_imag
            yr(k) = yr(k) + tmp_real
            yi(k) = yi(k) + tmp_imag
         next j
         k = k + n2
      loop until k >= n_points
      n1 = n1 - 1
      n2 = n2 / 2
   next i
! now do the bit reversal stuff
   for i = 1 to n_points
      p = p_val(i-1,n_power) + 1
      if p > i then     ! switch places
         tmp_real = yr(i)
         tmp_imag = yi(i)
         yr(i) = yr(p)
         yi(i) = yi(p)
         yr(p) = tmp_real
         yi(p) = tmp_imag
      end if
   next i
! finally have to repack results into the array y()
   for i = 1 to n_points
      if f_type > 0 then   ! a forward transform
         y(2*i-1) = yr(i)
         y(2*i) = yi(i)
      else                 ! a backward transform
         y(i) = 2 * yr(i) / n_points
      end if
   next i
end sub
! bit reversal function - does double duty
function p_val(a,b)
   tmp = int(a)
   p = 0
   for i = 1 to b
      p = 2 * p + tmp - 2 * int(tmp / 2)
      tmp = int(tmp / 2)
   next i
   p_val = p
end function