! this program uses the FFT to compute the power spectrum of a signal
! the data should be in a file as
! t1, y1
! t2, y2
! t3, y3
! etc.
! the power spectrum is written to a file
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program power
   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)               ! calculate corresponding frequencies
   call power(t,y)                          ! compute power spectrum
   call display(t,y)                        ! display and store results
   call save(t,y,store$)
end
! compute the power spectrum from the Fourier components in y()
sub power(t(),y())
   n = size(t) / 4
   for i = 1 to n
      y(i) = y(2*i-1)^2 + y(2*i)^2
   next i
   mat redim t(n),y(n)
end sub
! 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)
   dt = t(3) - t(1)      ! transforms slightly differently
   for i = 1 to n-1 
      t(i) = (i-1) / (n * dt)
   next i
end sub
! save the power spectrum
sub save(t(),y(),store$)
   open #1: name store$ & ".power", organization text, create new
   for i = 1 to size(t)
      print #1: t(i),",",y(i)  
   next i
   close #1
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$
   isign = 1
   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())
   call datagraph(t,d,3,1,"black")    ! graph power spectrum
   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