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;

}

@binkleym - Thank you for writing this, it was incredibly useful to me.
Also thank you @mbijon for making this code available!

Please let me know if either of you have a page where I can donate you some coffee/beer money.
Thanks!

In above Code What is an input ?

Hi @vasirajan. The code I wrote computes the FFT (Fast Fourier Transform) of a vector whose length is a power of 2 (ie, a vector whose length can be written as 2^n, where n is an integer). The algorithm is a straight port of the FFT code from "Numerical Algorithms in C".

You should be able to take the fft function above and do something like:

$array = array(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16);
$FFT_array = Fourier($array, 1); # Compute FFT of the array
$Inverse_FFT_array = Fourier($FFT_array, -1) # Compute inverse FFT, which should equal our original input vector.

What the input array is depends on what you're trying to do. FFT's are often used in signal analysis, image analysis, and several other fields.

And thanks @stampycode! Just do something nice for someone and tell the universe I did it. ;-)

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.