Fast Fourier Transform in PHP

fft.php

<?php

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

###################################################################
# PHP_Fourier 0.03b
# Original Fortran source by Numerical Recipies
# PHP port by Mathew Binkley ([email protected])
###################################################################

###################################################################
# Fourier($data, $sign) - Performs the FFT on the *complex* 
#                         array $data
#
# Presumes that count($data) is an integer power of two (ie: 2^n)
# (hint: When your $data length is not a power of 2, pad with zeros to the next-higher power.)
#
# $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:
#      FFT operates on an array 
#      $data = array(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16); # 16 = 2^4
#
#      Compute FFT of the `$data` array:
#      $FFT_array = Fourier($data, 1);
#
#      Compute inverse FFT, which should equal our original `$data` vector:
#      $Inverse_FFT_array = Fourier($FFT_array, -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;

}

You're welcome @stampycode, but all the work represented here is from @binkleym. I'm a big fan of buying a sandwich for someone down on their luck and hungry. Or whatever other nice thing you can think of to do.

I'm going to update the doc header in the Gist now. It will note the required n^2 array-length on the input. That made sense to me back in 2011 when I was using FFT regularly. But I would have needed time to remember that if not for Mathew's comment today.

@mbijon @binkleym

It looks like this function is not working properly.

I've skipped the normalization to compare results with NumPy FFT
//$data[$i] *= sqrt(2/$n);

for example:

$data = [4, 2, 5, 7, 6, 8, 9, 1];
print_r(Fourier($data, 1));

Array
(
    [0] => 24
    [1] => 18
    [2] => -8
    [3] => -10
    [4] => -4
    [5] => 2
    [6] => 4
    [7] => -2
)

and numpy:

import numpy as np
data = np.array([4 + 2j, 5 + 7j, 6 + 8j, 9 + 1j])
print(np.fft.fft(data))

[24.+18.j 4. -2.j -4. +2.j -8.-10.j]

All result values (except those in the beginning and in the middle) are wrong.
Only if there are 4 or less values in the array (2 or less complex numbers) there are no errors.