fft.php

<?php

// !!! Warning: for reference, not debugged

###################################################################
# PHP_Fourier 0.03
# Original Fortran source by Numerical Recipies
# PHP port by Mathew Binkley (binkleym@nukote.com)
###################################################################

###################################################################
# Fourier($data, $sign) - Performs the FFT on the *complex* 
#                         array $data
#
# Presumes that count($data) is an integer power of two
#
# $data[even] holds the real portion
# $data[odd] hold the imaginary portion
# 
# Example: (1 + 2i) ->  $data[0] = 1; $data[1] = 2;
#
# $sign = 1  performs the Fourier Transform
# $sign = -1 performs the Inverse Fourier Transform
#
# Use:
#      $fourier_array = Fourier($inputarray, 1);
#
###################################################################

function Fourier($input, $isign) {

   #####################################################################
   # We need to shift the array up one because this script is a direct 
   # port of the fortran program from NR.  Should fix in future.
   #####################################################################
   $data[0] = 0; 
   for ($i = 0; $i < count($input); $i++) $data[($i+1)] = $input[$i];

   $n = count($input);

   $j = 1;

   for ($i = 1; $i < $n; $i += 2) {
      if ($j > $i) {
         list($data[($j+0)], $data[($i+0)]) = array($data[($i+0)], $data[($j+0)]);
         list($data[($j+1)], $data[($i+1)]) = array($data[($i+1)], $data[($j+1)]);
      }

      $m = $n >> 1;

      while (($m >= 2) && ($j > $m)) {
         $j -= $m;
         $m = $m >> 1;
      }

      $j += $m;

   }

   $mmax = 2;
 
   while ($n > $mmax) {  # Outer loop executed log2(nn) times
      $istep = $mmax << 1;

      $theta = $isign * 2*pi()/$mmax;

      $wtemp = sin(0.5 * $theta);
      $wpr   = -2.0*$wtemp*$wtemp;
      $wpi   = sin($theta);
 
      $wr = 1.0;
      $wi = 0.0;
      for ($m = 1; $m < $mmax; $m += 2) {  # Here are the two nested inner loops
         for ($i = $m; $i <= $n; $i+= $istep) {

            $j = $i + $mmax;

            $tempr = $wr * $data[$j]     - $wi * $data[($j+1)];
            $tempi = $wr * $data[($j+1)] + $wi * $data[$j];

            $data[$j]     = $data[$i]     - $tempr;
            $data[($j+1)] = $data[($i+1)] - $tempi;

            $data[$i]     += $tempr;
            $data[($i+1)] += $tempi;

         }
         $wtemp = $wr;
         $wr = ($wr * $wpr) - ($wi    * $wpi) + $wr;
         $wi = ($wi * $wpr) + ($wtemp * $wpi) + $wi;
      }
      $mmax = $istep;
   }

   for ($i = 1; $i < count($data); $i++) { 
      $data[$i] *= sqrt(2/$n);                   # Normalize the data
      if (abs($data[$i]) < 1E-8) $data[$i] = 0;  # Let's round small numbers to zero
      $input[($i-1)] = $data[$i];                # We need to shift array back (see beginning)
   }

   return $input;

}

?>