-
-
Save BrianHannay/c5c23bcc2d471e40d4d6f22c85adf60b to your computer and use it in GitHub Desktop.
Fast Fourier Transform in PHP
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
<?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; | |
} | |
?> |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment