Skip to content

Instantly share code, notes, and snippets.

@mbijon
Last active February 17, 2024 03:46
  • Star 18 You must be signed in to star a gist
  • Fork 3 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save mbijon/1332348 to your computer and use it in GitHub Desktop.
Fast Fourier Transform in PHP
<?php
// !!! Warning: for reference, not debugged
###################################################################
# PHP_Fourier 0.03b
# 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 (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;
}
@vincentcox
Copy link

Hi there,
I'm really desperate at the moment. I have to make a FFT analyser on my raspberry pi to get a graphical representation of the sound it records. I spend the last 2 days messing with C++ libraries, but I can't compile it on my raspberry pi or it just dont work.
PhP is something I am more familiar with, and I stumbled over this code of you.
How do you use this piece of code?
Do I have to pass it an array with amplitues of a wave signal? What does it output?

@binkleym
Copy link

Hi Vincentcox. I'm the author of the original PHP FFT code above. The above code was written to "scratch an itch". I was trying to do FFT's to improve forecasting of some sales data.

If you're doing FFT's on a RasPi, I'd recommend starting here:

https://www.raspberrypi.org/blog/accelerating-fourier-transforms-using-the-gpu/

@vasirajan
Copy link

In above Code What is an input ?

@stampycode
Copy link

@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!

@binkleym
Copy link

binkleym commented Aug 1, 2019

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. ;-)

@mbijon
Copy link
Author

mbijon commented Aug 1, 2019

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.

@rst59
Copy link

rst59 commented Apr 4, 2023

@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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment