Hey, it's the holidays, and some of our mods won't be available as much as they'd like to be. Those of you with mod powers, now would be a good time to help step in and stem the constant flood of spam... --Michael Mol

Fast Fourier transform

From Rosetta Code
(Redirected from FFT)
Jump to: navigation, search
Task
Fast Fourier transform
You are encouraged to solve this task according to the task description, using any language you may know.

The purpose of this task is to calculate the FFT (Fast Fourier Transform) of an input sequence.
The most general case allows for complex numbers at the input and results in a sequence of equal length, again of complex numbers. If you need to restrict yourself to real numbers, the output should be the magnitude (i.e. sqrt(re²+im²)) of the complex result.

The classic version is the recursive Cooley–Tukey FFT. Wikipedia has pseudocode for that. Further optimizations are possible but not required.

Contents

[edit] Ada

The FFT function is defined as a generic function, instantiated upon a user instance of Ada.Numerics.Generic_Complex_Arrays.

 
with Ada.Numerics.Generic_Complex_Arrays;
 
generic
with package Complex_Arrays is
new Ada.Numerics.Generic_Complex_Arrays (<>);
use Complex_Arrays;
function Generic_FFT (X : Complex_Vector) return Complex_Vector;
 
 
with Ada.Numerics;
with Ada.Numerics.Generic_Complex_Elementary_Functions;
 
function Generic_FFT (X : Complex_Vector) return Complex_Vector is
 
package Complex_Elementary_Functions is
new Ada.Numerics.Generic_Complex_Elementary_Functions
(Complex_Arrays.Complex_Types);
 
use Ada.Numerics;
use Complex_Elementary_Functions;
use Complex_Arrays.Complex_Types;
 
function FFT (X : Complex_Vector; N, S : Positive)
return Complex_Vector is
begin
if N = 1 then
return (1..1 => X (X'First));
else
declare
F : constant Complex  := exp (Pi * j / Real_Arrays.Real (N/2));
Even : Complex_Vector := FFT (X, N/2, 2*S);
Odd  : Complex_Vector := FFT (X (X'First + S..X'Last), N/2, 2*S);
begin
for K in 0..N/2 - 1 loop
declare
T : constant Complex := Odd (Odd'First + K) / F ** K;
begin
Odd (Odd'First + K) := Even (Even'First + K) - T;
Even (Even'First + K) := Even (Even'First + K) + T;
end;
end loop;
return Even & Odd;
end;
end if;
end FFT;
begin
return FFT (X, X'Length, 1);
end Generic_FFT;
 

Example:

 
with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays;
with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;
with Ada.Text_IO; use Ada.Text_IO;
 
with Ada.Numerics.Complex_Elementary_Functions;
with Generic_FFT;
 
procedure Example is
function FFT is new Generic_FFT (Ada.Numerics.Complex_Arrays);
X : Complex_Vector := (1..4 => (1.0, 0.0), 5..8 => (0.0, 0.0));
Y : Complex_Vector := FFT (X);
begin
Put_Line (" X FFT X ");
for I in Y'Range loop
Put (X (I - Y'First + X'First), Aft => 3, Exp => 0);
Put (" ");
Put (Y (I), Aft => 3, Exp => 0);
New_Line;
end loop;
end;
 
Output:
       X              FFT X 
( 1.000, 0.000) ( 4.000, 0.000)
( 1.000, 0.000) ( 1.000,-2.414)
( 1.000, 0.000) ( 0.000, 0.000)
( 1.000, 0.000) ( 1.000,-0.414)
( 0.000, 0.000) ( 0.000, 0.000)
( 0.000, 0.000) ( 1.000, 0.414)
( 0.000, 0.000) ( 0.000, 0.000)
( 0.000, 0.000) ( 1.000, 2.414)

[edit] ALGOL 68

Translation of: Python
Note: This specimen retains the original Python coding style.
Works with: ALGOL 68 version Revision 1 - one minor extension to language used - PRAGMA READ, similar to C's #include directive.
Works with: ALGOL 68G version Any - tested with release algol68g-2.3.5.
File: Template.Fast_Fourier_transform.a68
PRIO DICE = 9; # ideally = 11 #
 
OP DICE = ([]SCALAR in, INT step)[]SCALAR: (
### Dice the array, extract array values a "step" apart ###
IF step = 1 THEN
in
ELSE
INT upb out := 0;
[(UPB in-LWB in)%step+1]SCALAR out;
FOR index FROM LWB in BY step TO UPB in DO
out[upb out+:=1] := in[index] OD;
out[@LWB in]
FI
);
 
PROC fft = ([]SCALAR in t)[]SCALAR: (
### The Cooley-Tukey FFT algorithm ###
IF LWB in t >= UPB in t THEN
in t[@0]
ELSE
[]SCALAR t = in t[@0];
INT n = UPB t + 1, half n = n % 2;
[LWB t:UPB t]SCALAR coef;
 
[]SCALAR even = fft(t DICE 2),
odd = fft(t[1:]DICE 2);
 
COMPL i = 0 I 1;
 
REAL w = 2*pi / n;
FOR k FROM LWB t TO half n-1 DO
COMPL cis t = scalar exp(0 I (-w * k))*odd[k];
coef[k] := even[k] + cis t;
coef[k + half n] := even[k] - cis t
OD;
coef
FI
);
File: test.Fast_Fourier_transform.a68
#!/usr/local/bin/a68g --script #
# -*- coding: utf-8 -*- #
 
MODE SCALAR = COMPL;
PROC (COMPL)COMPL scalar exp = complex exp;
PR READ "Template.Fast_Fourier_transform.a68" PR
 
FORMAT real fmt := $g(0,3)$;
FORMAT real array fmt := $f(real fmt)", "$;
FORMAT compl fmt := $f(real fmt)"⊥"f(real fmt)$;
FORMAT compl array fmt := $f(compl fmt)", "$;
 
test:(
[]COMPL
tooth wave ft = fft((1, 1, 1, 1, 0, 0, 0, 0)),
one and a quarter wave ft = fft((0, 0.924, 0.707,-0.383,-1,-0.383, 0.707, 0.924,
0,-0.924,-0.707, 0.383, 1, 0.383,-0.707,-0.924));
printf((
$"Tooth wave: "$,compl array fmt, tooth wave ft, $l$,
$"1¼ cycle wave: "$, compl array fmt, one and a quarter wave ft, $l$
))
)
Output:
Tooth wave: 4.000⊥.000, 1.000⊥-2.414, .000⊥.000, 1.000⊥-.414, .000⊥.000, 1.000⊥.414, .000⊥.000, 1.000⊥2.414, 
1¼ cycle wave: .000⊥.000, .000⊥.001, .000⊥.000, .000⊥-8.001, .000⊥.000, -.000⊥-.001, .000⊥.000, .000⊥.001, .000⊥.000, .000⊥-.001, .000⊥.000, -.000⊥.001, .000⊥.000, -.000⊥8.001, .000⊥.000, -.000⊥-.001, 

[edit] BBC BASIC

      @% = &60A
 
DIM Complex{r#, i#}
DIM in{(7)} = Complex{}, out{(7)} = Complex{}
DATA 1, 1, 1, 1, 0, 0, 0, 0
 
PRINT "Input (real, imag):"
FOR I% = 0 TO 7
READ in{(I%)}.r#
PRINT in{(I%)}.r# "," in{(I%)}.i#
NEXT
 
PROCfft(out{()}, in{()}, 0, 1, DIM(in{()},1)+1)
 
PRINT "Output (real, imag):"
FOR I% = 0 TO 7
PRINT out{(I%)}.r# "," out{(I%)}.i#
NEXT
END
 
DEF PROCfft(b{()}, o{()}, B%, S%, N%)
LOCAL I%, t{} : DIM t{} = Complex{}
IF S% < N% THEN
PROCfft(o{()}, b{()}, B%, S%*2, N%)
PROCfft(o{()}, b{()}, B%+S%, S%*2, N%)
FOR I% = 0 TO N%-1 STEP 2*S%
t.r# = COS(-PI*I%/N%)
t.i# = SIN(-PI*I%/N%)
PROCcmul(t{}, o{(B%+I%+S%)})
b{(B%+I% DIV 2)}.r# = o{(B%+I%)}.r# + t.r#
b{(B%+I% DIV 2)}.i# = o{(B%+I%)}.i# + t.i#
b{(B%+(I%+N%) DIV 2)}.r# = o{(B%+I%)}.r# - t.r#
b{(B%+(I%+N%) DIV 2)}.i# = o{(B%+I%)}.i# - t.i#
NEXT
ENDIF
ENDPROC
 
DEF PROCcmul(c{},d{})
LOCAL r#, i#
r# = c.r#*d.r# - c.i#*d.i#
i# = c.r#*d.i# + c.i#*d.r#
c.r# = r#
c.i# = i#
ENDPROC
 
Output:
Input (real, imag):
         1,         0
         1,         0
         1,         0
         1,         0
         0,         0
         0,         0
         0,         0
         0,         0
Output (real, imag):
         4,         0
         1,  -2.41421
         0,         0
         1, -0.414214
         0,         0
         1,  0.414214
         0,         0
         1,   2.41421

[edit] C

Inplace FFT with O(n) memory usage. Note: array size is assumed to be power of 2 and not checked by code; you can just pad it with 0 otherwise.

Also, complex is C99 standard.
 
 
#include <stdio.h>
#include <math.h>
#include <complex.h>
 
double PI;
typedef double complex cplx;
 
void _fft(cplx buf[], cplx out[], int n, int step)
{
if (step < n) {
_fft(out, buf, n, step * 2);
_fft(out + step, buf + step, n, step * 2);
 
for (int i = 0; i < n; i += 2 * step) {
cplx t = cexp(-I * PI * i / n) * out[i + step];
buf[i / 2] = out[i] + t;
buf[(i + n)/2] = out[i] - t;
}
}
}
 
void fft(cplx buf[], int n)
{
cplx out[n];
for (int i = 0; i < n; i++) out[i] = buf[i];
 
_fft(buf, out, n, 1);
}
 
 
void show(const char * s, cplx buf[]) {
printf("%s", s);
for (int i = 0; i < 8; i++)
if (!cimag(buf[i]))
printf("%g ", creal(buf[i]));
else
printf("(%g, %g) ", creal(buf[i]), cimag(buf[i]));
}
 
int main()
{
PI = atan2(1, 1) * 4;
cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};
 
show("Data: ", buf);
fft(buf, 8);
show("\nFFT : ", buf);
 
return 0;
}
 
 
Output:
Data: 1 1 1 1 0 0 0 0 
FFT : 4 (1, -2.41421) 0 (1, -0.414214) 0 (1, 0.414214) 0 (1, 2.41421)

[edit] C++

#include <complex>
#include <iostream>
#include <valarray>
 
const double PI = 3.141592653589793238460;
 
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;
 
// Cooley–Tukey FFT (in-place)
void fft(CArray& x)
{
const size_t N = x.size();
if (N <= 1) return;
 
// divide
CArray even = x[std::slice(0, N/2, 2)];
CArray odd = x[std::slice(1, N/2, 2)];
 
// conquer
fft(even);
fft(odd);
 
// combine
for (size_t k = 0; k < N/2; ++k)
{
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
 
// inverse fft (in-place)
void ifft(CArray& x)
{
// conjugate the complex numbers
x = x.apply(std::conj);
 
// forward fft
fft( x );
 
// conjugate the complex numbers again
x = x.apply(std::conj);
 
// scale the numbers
x /= x.size();
}
 
int main()
{
const Complex test[] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
CArray data(test, 8);
 
// forward fft
fft(data);
 
std::cout << "fft" << std::endl;
for (int i = 0; i < 8; ++i)
{
std::cout << data[i] << std::endl;
}
 
// inverse fft
ifft(data);
 
std::cout << std::endl << "ifft" << std::endl;
for (int i = 0; i < 8; ++i)
{
std::cout << data[i] << std::endl;
}
return 0;
}
Output:
fft
(4,0)
(1,-2.41421)
(0,0)
(1,-0.414214)
(0,0)
(1,0.414214)
(0,0)
(1,2.41421)

ifft
(1,-0)
(1,-5.55112e-017)
(1,2.4895e-017)
(1,-5.55112e-017)
(5.55112e-017,0)
(5.55112e-017,5.55112e-017)
(0,-2.4895e-017)
(5.55112e-017,5.55112e-017)

[edit] D

[edit] Standard Version

void main() {
import std.stdio, std.numeric;
 
[1.0, 1, 1, 1, 0, 0, 0, 0].fft.writeln;
}
Output:
[4+0i, 1-2.41421i, 0-0i, 1-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

[edit] creals Version

Built-in complex numbers will be deprecated.

import std.stdio, std.algorithm, std.range, std.math;
 
const(creal)[] fft(in creal[] x) /*pure nothrow*/ {
immutable N = x.length;
if (N <= 1) return x;
const ev = x.stride(2).array.fft;
const od = x[1 .. $].stride(2).array.fft;
auto l = iota(N / 2).map!(k => ev[k] + expi(-2*PI * k/N) * od[k]);
auto r = iota(N / 2).map!(k => ev[k] - expi(-2*PI * k/N) * od[k]);
return l.chain(r).array;
}
 
void main() {
[1.0L+0i, 1, 1, 1, 0, 0, 0, 0].fft.writeln;
}
Output:
[4+0i, 1+-2.41421i, 0+0i, 1+-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

[edit] Phobos Complex Version

import std.stdio, std.algorithm, std.range, std.math, std.complex;
 
auto fft(T)(in T[] x) /*pure nothrow*/ {
immutable N = x.length;
if (N <= 1) return x;
const ev = x.stride(2).array.fft;
const od = x[1 .. $].stride(2).array.fft;
alias E = std.complex.expi;
auto l = iota(N / 2).map!(k => ev[k] + T(E(-2* PI * k/N)) * od[k]);
auto r = iota(N / 2).map!(k => ev[k] - T(E(-2* PI * k/N)) * od[k]);
return l.chain(r).array;
}
 
void main() {
[1.0, 1, 1, 1, 0, 0, 0, 0].map!complex.array.fft.writeln;
}
Output:
[4+0i, 1-2.41421i, 0+0i, 1-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

[edit] Factor

 
IN: USE math.transforms.fft
IN: { 1 1 1 1 0 0 0 0 } fft .
{
C{ 4.0 0.0 }
C{ 1.0 -2.414213562373095 }
C{ 0.0 0.0 }
C{ 1.0 -0.4142135623730949 }
C{ 0.0 0.0 }
C{ 0.9999999999999999 0.4142135623730949 }
C{ 0.0 0.0 }
C{ 0.9999999999999997 2.414213562373095 }
}
 

[edit] Fortran

 
module fft_mod
implicit none
integer, parameter :: dp=selected_real_kind(15,300)
real(kind=dp), parameter :: pi=3.141592653589793238460_dp
contains
 
! In place Cooley-Tukey FFT
recursive subroutine fft(x)
complex(kind=dp), dimension(:), intent(inout) :: x
complex(kind=dp) :: t
integer :: N
integer :: i
complex(kind=dp), dimension(:), allocatable :: even, odd
 
N=size(x)
 
if(N .le. 1) return
 
allocate(odd((N+1)/2))
allocate(even(N/2))
 
! divide
odd =x(1:N:2)
even=x(2:N:2)
 
! conquer
call fft(odd)
call fft(even)
 
! combine
do i=1,N/2
 
 
t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),KIND=DP))*even(i)
x(i) = odd(i) + t
x(i+N/2) = odd(i) - t
end do
 
deallocate(odd)
deallocate(even)
 
end subroutine fft
 
end module fft_mod
 
program test
use fft_mod
implicit none
complex(kind=dp), dimension(8) :: data = (/1.0, 1.0, 1.0, 1.0, 0.0,
 
0.0, 0.0, 0.0/)
integer :: i
 
call fft(data)
 
do i=1,8
write(*,'("(", F20.15, ",", F20.15, "i )")') data(i)
end do
 
end program test
Output:
(   4.000000000000000,   0.000000000000000i )
(   1.000000000000000,  -2.414213562373095i )
(   0.000000000000000,   0.000000000000000i )
(   1.000000000000000,  -0.414213562373095i )
(   0.000000000000000,   0.000000000000000i )
(   1.000000000000000,   0.414213562373095i )
(   0.000000000000000,   0.000000000000000i )
(   1.000000000000000,   2.414213562373095i )

[edit] GAP

# Here an implementation with no optimization (O(n^2)).
# In GAP, E(n) = exp(2*i*pi/n), a primitive root of the unity.
 
Fourier := function(a)
local n, z;
n := Size(a);
z := E(n);
return List([0 .. n - 1], k -> Sum([0 .. n - 1], j -> a[j + 1]*z^(-k*j)));
end;
 
InverseFourier := function(a)
local n, z;
n := Size(a);
z := E(n);
return List([0 .. n - 1], k -> Sum([0 .. n - 1], j -> a[j + 1]*z^(k*j)))/n;
end;
 
Fourier([1, 1, 1, 1, 0, 0, 0, 0]);
# [ 4, 1-E(8)-E(8)^2-E(8)^3, 0, 1-E(8)+E(8)^2-E(8)^3,
# 0, 1+E(8)-E(8)^2+E(8)^3, 0, 1+E(8)+E(8)^2+E(8)^3 ]
 
InverseFourier(last);
# [ 1, 1, 1, 1, 0, 0, 0, 0 ]

[edit] Go

package main
 
import (
"fmt"
"math"
"math/cmplx"
)
 
func ditfft2(x []float64, y []complex128, n, s int) {
if n == 1 {
y[0] = complex(x[0], 0)
return
}
ditfft2(x, y, n/2, 2*s)
ditfft2(x[s:], y[n/2:], n/2, 2*s)
for k := 0; k < n/2; k++ {
tf := cmplx.Rect(1, -2*math.Pi*float64(k)/float64(n)) * y[k+n/2]
y[k], y[k+n/2] = y[k]+tf, y[k]-tf
}
}
 
func main() {
x := []float64{1, 1, 1, 1, 0, 0, 0, 0}
y := make([]complex128, len(x))
ditfft2(x, y, len(x), 1)
for _, c := range y {
fmt.Printf("%8.4f\n", c)
}
}
Output:
(  4.0000 +0.0000i)
(  1.0000 -2.4142i)
(  0.0000 +0.0000i)
(  1.0000 -0.4142i)
(  0.0000 +0.0000i)
(  1.0000 +0.4142i)
(  0.0000 +0.0000i)
(  1.0000 +2.4142i)

[edit] Haskell

import Data.Complex
 
-- Cooley-Tukey
fft [] = []
fft [x] = [x]
fft xs = zipWith (+) ys ts ++ zipWith (-) ys ts
where n = length xs
ys = fft evens
zs = fft odds
(evens, odds) = split xs
split [] = ([], [])
split [x] = ([x], [])
split (x:y:xs) = (x:xt, y:yt) where (xt, yt) = split xs
ts = zipWith (\z k -> exp' k n * z) zs [0..]
exp'
k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)
 
main = mapM_ print $ fft [1,1,1,1,0,0,0,0]
Output:
4.0 :+ 0.0
1.0 :+ (-2.414213562373095)
0.0 :+ 0.0
1.0 :+ (-0.4142135623730949)
0.0 :+ 0.0
0.9999999999999999 :+ 0.4142135623730949
0.0 :+ 0.0
0.9999999999999997 :+ 2.414213562373095

[edit] J

Based on j:Essays/FFT, with some simplifications -- sacrificing accuracy, optimizations and convenience which are not relevant to the task requirements, for clarity:

cube  =: ($~ q:@#) :. ,
rou =: ^@j.@o.@(% #)@i.@-: NB. roots of unity
floop =: 4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.'
fft =: ] floop&.cube rou@#

Example (first row of result is sine, second row of result is fft of the first row, (**+)&.+. cleans an irrelevant least significant bit of precision from the result so that it displays nicely):

   (**+)&.+. (,: fft) 1 o. 2p1*3r16 * i.16
0 0.92388 0.707107 0.382683 1 0.382683 0.707107 0.92388 0 0.92388 0.707107 0.382683 1 0.382683 0.707107 0.92388
0 0 0 0j8 0 0 0 0 0 0 0 0 0 0j8 0 0

Here is a representation of an example which appears in some of the other implementations, here:

   Re=: {.@+.@fft
Im=: {:@+.@fft
M=: 4#1 0
M
1 1 1 1 0 0 0 0
Re M
4 1 0 1 0 1 0 1
Im M
0 2.41421 0 0.414214 0 _0.414214 0 _2.41421

Note that Re and Im are not functions of 1 and 0 but are functions of the complete sequence.

Also note that J uses a different character for negative sign than for subtraction, to eliminate ambiguity (is this a single list of numbers or are lists being subtracted?).

[edit] JavaScript

Complex fourier transform & it's inverse reimplemented from the C++ & Python variants on this page.

/*
complex fast fourier transform and inverse from
http://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B
*/

function icfft(amplitudes)
{
var N = amplitudes.length;
var iN = 1 / N;
 
//conjugate if imaginary part is not 0
for(var i = 0 ; i < N; ++i)
if(amplitudes[i] instanceof Complex)
amplitudes[i].im = -amplitudes[i].im;
 
//apply fourier transform
amplitudes = cfft(amplitudes)
 
for(var i = 0 ; i < N; ++i)
{
//conjugate again
amplitudes[i].im = -amplitudes[i].im;
//scale
amplitudes[i].re *= iN;
amplitudes[i].im *= iN;
}
return amplitudes;
}
 
function cfft(amplitudes)
{
var N = amplitudes.length;
if( N <= 1 )
return amplitudes;
 
var hN = N / 2;
var even = [];
var odd = [];
even.length = hN;
odd.length = hN;
for(var i = 0; i < hN; ++i)
{
even[i] = amplitudes[i*2];
odd[i] = amplitudes[i*2+1];
}
even = cfft(even);
odd = cfft(odd);
 
var a = -2*Math.PI;
for(var k = 0; k < hN; ++k)
{
if(!(even[k] instanceof Complex))
even[k] = new Complex(even[k], 0);
if(!(odd[k] instanceof Complex))
odd[k] = new Complex(odd[k], 0);
var p = k/N;
var t = new Complex(0, a * p);
t.cexp(t).mul(odd[k], t);
amplitudes[k] = even[k].add(t, odd[k]);
amplitudes[k + hN] = even[k].sub(t, even[k]);
}
return amplitudes;
}
 
//test code
//console.log( cfft([1,1,1,1,0,0,0,0]) );
//console.log( icfft(cfft([1,1,1,1,0,0,0,0])) );

Very very basic Complex number that provides only the components required by the code above.

/*
basic complex number arithmetic from
http://rosettacode.org/wiki/Fast_Fourier_transform#Scala
*/

function Complex(re, im)
{
this.re = re;
this.im = im || 0.0;
}
Complex.prototype.add = function(other, dst)
{
dst.re = this.re + other.re;
dst.im = this.im + other.im;
return dst;
}
Complex.prototype.sub = function(other, dst)
{
dst.re = this.re - other.re;
dst.im = this.im - other.im;
return dst;
}
Complex.prototype.mul = function(other, dst)
{
//cache re in case dst === this
var r = this.re * other.re - this.im * other.im;
dst.im = this.re * other.im + this.im * other.re;
dst.re = r;
return dst;
}
Complex.prototype.cexp = function(dst)
{
var er = Math.exp(this.re);
dst.re = er * Math.cos(this.im);
dst.im = er * Math.sin(this.im);
return dst;
}
Complex.prototype.log = function()
{
/*
although 'It's just a matter of separating out the real and imaginary parts of jw.' is not a helpful quote
the actual formula I found here and the rest was just fiddling / testing and comparing with correct results.
http://cboard.cprogramming.com/c-programming/89116-how-implement-complex-exponential-functions-c.html#post637921
*/

if( !this.re )
console.log(this.im.toString()+'j');
else if( this.im < 0 )
console.log(this.re.toString()+this.im.toString()+'j');
else
console.log(this.re.toString()+'+'+this.im.toString()+'j');
}

[edit] jq

Currently jq has no support for complex numbers, so the following implementation uses [x,y] to represent the complex number x+iy.

[edit] Complex number arithmetic

 
# multiplication of real or complex numbers
def cmult(x; y):
if (x|type) == "number" then
if (y|type) == "number" then [ x*y, 0 ]
else [x * y[0], x * y[1]]
end
elif (y|type) == "number" then cmult(y;x)
else [ x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]]
end;
 
def cplus(x; y):
if (x|type) == "number" then
if (y|type) == "number" then [ x+y, 0 ]
else [ x + y[0], y[1]]
end
elif (y|type) == "number" then cplus(y;x)
else [ x[0] + y[0], x[1] + y[1] ]
end;
 
def cminus(x; y): cplus(x; cmult(-1; y));
 
# e(ix) = cos(x) + i sin(x)
def expi(x): [ (x|cos), (x|sin) ];

[edit] FFT

def fft:
length as $N
| if $N <= 1 then .
else ( [ .[ range(0; $N; 2) ] ] | fft) as $even
| ( [ .[ range(1; $N; 2) ] ] | fft) as $odd
| (1|atan * 4) as $pi
| [ range(0; $N/2) | cplus($even[.]; cmult( expi(-2*$pi*./$N); $odd[.] )) ] +
[ range(0; $N/2) | cminus($even[.]; cmult( expi(-2*$pi*./$N); $odd[.] )) ]
end;

Example:

[1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0] | fft
 
Output:
[[4,-0],[1,-2.414213562373095],
 [0,0],[1,-0.4142135623730949],
 [0,0],[0.9999999999999999,0.4142135623730949],
 [0,0],[0.9999999999999997,2.414213562373095]]

[edit] Julia

Julia has a built-in FFT function:

fft([1,1,1,1,0,0,0,0])
Output:
8-element Array{Complex{Float64},1}:
4.0+0.0im
1.0-2.41421im
0.0+0.0im
1.0-0.414214im
0.0+0.0im
1.0+0.414214im
0.0+0.0im
1.0+2.41421im

[edit] Liberty BASIC

 
P =8
S =int( log( P) /log( 2) +0.9999)
 
Pi =3.14159265
R1 =2^S
 
R =R1 -1
R2 =div( R1, 2)
R4 =div( R1, 4)
R3 =R4 +R2
 
Dim Re( R1), Im( R1), Co( R3)
 
for N =0 to P -1
read dummy: Re( N) =dummy
read dummy: Im( N) =dummy
next N
 
data 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0
 
S2 =div( S, 2)
S1 =S -S2
P1 =2^S1
P2 =2^S2
 
dim V( P1 -1)
V( 0) =0
DV =1
DP =P1
 
for J =1 to S1
HA =div( DP, 2)
PT =P1 -HA
for I =HA to PT step DP
V( I) =V( I -HA) +DV
next I
DV =DV +DV
DP =HA
next J
 
K =2 *Pi /R1
 
for X =0 to R4
COX =cos( K *X)
Co( X) =COX
Co( R2 -X) =0 -COX
Co( R2 +X) =0 -COX
next X
 
print "FFT: bit reversal"
 
for I =0 to P1 -1
IP =I *P2
for J =0 to P2 -1
H =IP +J
G =V( J) *P2 +V( I)
if G >H then temp =Re( G): Re( G) =Re( H): Re( H) =temp
if G >H then temp =Im( G): Im( G) =Im( H): Im( H) =temp
next J
next I
 
T =1
 
for stage =0 to S -1
print " Stage:- "; stage
D =div( R2, T)
for Z =0 to T -1
L =D *Z
LS =L +R4
for I =0 to D -1
A =2 *I *T +Z
B =A +T
F1 =Re( A)
F2 =Im( A)
P1 =Co( L) *Re( B)
P2 =Co( LS) *Im( B)
P3 =Co( LS) *Re( B)
P4 =Co( L) *Im( B)
Re( A) =F1 +P1 -P2
Im( A) =F2 +P3 +P4
Re( B) =F1 -P1 +P2
Im( B) =F2 -P3 -P4
next I
next Z
T =T +T
next stage
 
print " M Re( M) Im( M)"
 
for M =0 to R
if abs( Re( M)) <10^-5 then Re( M) =0
if abs( Im( M)) <10^-5 then Im( M) =0
print " "; M, Re( M), Im( M)
next M
 
end
 
 
wait
 
function div( a, b)
div =int( a /b)
end function
 
end
 
   M          Re( M)       Im( M)
  0          4             0
  1          1.0           -2.41421356
  2          0             0
  3          1.0           -0.41421356
  4          0             0
  5          1.0           0.41421356
  6          0             0
  7          1.0           2.41421356

[edit] Maple

Maple has a built-in package DiscreteTransforms, and FourierTransform and InverseFourierTransform are in the commands available from that package. The FourierTransform command offers an FFT method by default.

 
with( DiscreteTransforms ):
 
FourierTransform( <1,1,1,1,0,0,0,0>, normalization=none );
 
                         [       4. + 0. I        ]
                         [                        ]
                         [1. - 2.41421356237309 I ]
                         [                        ]
                         [       0. + 0. I        ]
                         [                        ]
                         [1. - 0.414213562373095 I]
                         [                        ]
                         [       0. + 0. I        ]
                         [                        ]
                         [1. + 0.414213562373095 I]
                         [                        ]
                         [       0. + 0. I        ]
                         [                        ]
                         [1. + 2.41421356237309 I ]

Optionally, the FFT may be performed inplace on a Vector of hardware double-precision complex floats.

 
v := Vector( [1,1,1,1,0,0,0,0], datatype=complex[8] ):
 
FourierTransform( v, normalization=none, inplace ):
 
v;
 
                         [       4. + 0. I        ]
                         [                        ]
                         [1. - 2.41421356237309 I ]
                         [                        ]
                         [       0. + 0. I        ]
                         [                        ]
                         [1. - 0.414213562373095 I]
                         [                        ]
                         [       0. + 0. I        ]
                         [                        ]
                         [1. + 0.414213562373095 I]
                         [                        ]
                         [       0. + 0. I        ]
                         [                        ]
                         [1. + 2.41421356237309 I ]
 
InverseFourierTransform( v, normalization=full, inplace ):
 
v;
 
                       [          1. + 0. I          ]
                       [                             ]
                       [          1. + 0. I          ]
                       [                             ]
                       [          1. + 0. I          ]
                       [                             ]
                       [          1. + 0. I          ]
                       [                             ]
                       [          0. + 0. I          ]
                       [                             ]
                       [          0. + 0. I          ]
                       [                             ]
                       [                   -17       ]
                       [5.55111512312578 10    + 0. I]
                       [                             ]
                       [          0. + 0. I          ]

[edit] Mathematica

Mathematica has a built-in FFT function which uses a proprietary algorithm developed at Wolfram Research. It also has an option to tune the algorithm for specific applications. The options shown below, while not default, produce output that is consistent with most other FFT routines.

 
Fourier[{1,1,1,1,0,0,0,0}, FourierParameters->{1,-1}]
 
Output:
{4. + 0. I, 1. - 2.4142136 I, 0. + 0. I, 1. - 0.41421356 I, 0. + 0. I, 1. + 0.41421356 I, 0. + 0. I, 1. + 2.4142136 I}

[edit] MATLAB / Octave

Matlab/Octave have a builtin FFT function.

 fft([1,1,1,1,0,0,0,0]')
 
Output:
ans =

   4.00000 + 0.00000i
   1.00000 - 2.41421i
   0.00000 + 0.00000i
   1.00000 - 0.41421i
   0.00000 + 0.00000i
   1.00000 + 0.41421i
   0.00000 - 0.00000i
   1.00000 + 2.41421i

[edit] Maxima

load(fft)$
fft([1, 2, 3, 4]);
[2.5, -0.5 * %i - 0.5, -0.5, 0.5 * %i - 0.5]

[edit] Nimrod

Translation of: Python
import math, complex, strutils
 
proc toComplex(x: float): TComplex = result.re = x
proc toComplex(x: TComplex): TComplex = x
 
# Works with floats and complex numbers as input
proc fft[T](x: openarray[T]): seq[TComplex] =
let n = x.len
result = newSeq[TComplex]()
if n <= 1:
for v in x: result.add toComplex(v)
return
var evens, odds = newSeq[T]()
for i, v in x:
if i mod 2 == 0: evens.add v
else: odds.add v
var (even, odd) = (fft(evens), fft(odds))
 
for k in 0 .. < n div 2:
result.add(even[k] + exp((0.0, -2*pi*float(k)/float(n))) * odd[k])
 
for k in 0 .. < n div 2:
result.add(even[k] - exp((0.0, -2*pi*float(k)/float(n))) * odd[k])
 
for i in fft(@[1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]):
echo formatFloat(abs(i), ffDecimal, 3)
Output:
4.000
2.613
-0.000
1.082
-0.000
1.082
-0.000
2.613

[edit] OCaml

This is a simple implementation of the Cooley-Tukey pseudo-code

open Complex
 
let fac k n =
let m2pi = -4.0 *. acos 0.0 in
polar 1.0 (m2pi*.(float k)/.(float n))
 
let merge l r n =
let f (k,t) x = (succ k, (mul (fac k n) x) :: t) in
let z = List.rev (snd (List.fold_left f (0,[]) r)) in
(List.map2 add l z) @ (List.map2 sub l z)
 
let fft lst =
let rec ditfft2 a n s =
if n = 1 then [List.nth lst a] else
let odd = ditfft2 a (n/2) (2*s) in
let even = ditfft2 (a+s) (n/2) (2*s) in
merge odd even n in
ditfft2 0 (List.length lst) 1;;
 
let show l =
let pr x = Printf.printf "(%f %f) " x.re x.im in
(List.iter pr l; print_newline ()) in
let indata = [one;one;one;one;zero;zero;zero;zero] in
show indata;
show (fft indata)
Output:
(1.000000 0.000000) (1.000000 0.000000) (1.000000 0.000000) (1.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) 
(4.000000 0.000000) (1.000000 -2.414214) (0.000000 0.000000) (1.000000 -0.414214) (0.000000 0.000000) (1.000000 0.414214) (0.000000 0.000000) (1.000000 2.414214) 

[edit] PARI/GP

Naive implementation, using the same testcase as Ada:

FFT(v)=my(t=-2*Pi*I/#v,tt);vector(#v,k,tt=t*(k-1);sum(n=0,#v-1,v[n+1]*exp(tt*n)));
FFT([1,1,1,1,0,0,0,0])
Output:
[4.0000000000000000000000000000000000000, 1.0000000000000000000000000000000000000 - 2.4142135623730950488016887242096980786*I, 0.E-37 + 0.E-38*I, 1.0000000000000000000000000000000000000 - 0.41421356237309504880168872420969807856*I, 0.E-38 + 0.E-37*I, 0.99999999999999999999999999999999999997 + 0.41421356237309504880168872420969807860*I, 4.701977403289150032 E-38 + 0.E-38*I, 0.99999999999999999999999999999999999991 + 2.4142135623730950488016887242096980785*I]

[edit] Perl

Translation of: Perl 6
use strict;
use warnings;
use Math::Complex;
 
sub fft {
return @_ if @_ == 1;
my @evn = fft(@_[grep { not $_ % 2 } 0 .. $#_ ]);
my @odd = fft(@_[grep { $_ % 2 } 1 .. $#_ ]);
my $twd = 2*i* pi / @_;
$odd[$_] *= exp( $_ * $twd ) for 0 .. $#odd;
return
(map { $evn[$_] + $odd[$_] } 0 .. $#evn ),
(map { $evn[$_] - $odd[$_] } 0 .. $#evn );
}
 
 
my @seq = 0 .. 15;
my $cycles = 3;
my @wave = map { sin( $_ * 2*pi/ @seq * $cycles ) } @seq;
print "wave: ", join " ", map { sprintf "%7.3f", $_ } @wave;
print "\n";
print "fft: ", join " ", map { sprintf "%7.3f", abs $_ } fft(@wave);
Output:
wave:   0.000   0.924   0.707  -0.383  -1.000  -0.383   0.707   0.924   0.000  -0.924  -0.707   0.383   1.000   0.383  -0.707  -0.924
fft:    0.000   0.000   0.000   8.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   8.000   0.000   0.000

[edit] Perl 6

sub fft {
return @_ if @_ == 1;
my @evn = fft( @_[0, 2 ...^ * >= @_] );
my @odd = fft( @_[1, 3 ...^ * >= @_] ) Z*
(1, * * cis( 2 * pi / @_ ) ... *);
return @evn »+« @odd, @evn »-« @odd;
}
 
my @seq = ^16;
my $cycles = 3;
my @wave = map { sin( 2*pi * $_ / @seq * $cycles ) }, @seq;
say "wave: ", @wave.fmt("%7.3f");
 
say "fft: ", fft(@wave)».abs.fmt("%7.3f");
Output:
wave:   0.000   0.924   0.707  -0.383  -1.000  -0.383   0.707   0.924   0.000  -0.924  -0.707   0.383   1.000   0.383  -0.707  -0.924
fft:    0.000   0.000   0.000   8.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   8.000   0.000   0.000

[edit] PicoLisp

Works with: PicoLisp version 3.1.0.3
# apt-get install libfftw3-dev
 
(scl 4)
 
(de FFTW_FORWARD . -1)
(de FFTW_ESTIMATE . 64)
 
(de fft (Lst)
(let
(Len (length Lst)
In (native "libfftw3.so" "fftw_malloc" 'N (* Len 16))
Out (native "libfftw3.so" "fftw_malloc" 'N (* Len 16))
P (native "libfftw3.so" "fftw_plan_dft_1d" 'N
Len In Out FFTW_FORWARD FFTW_ESTIMATE ) )
(struct In NIL (cons 1.0 (apply append Lst)))
(native "libfftw3.so" "fftw_execute" NIL P)
(prog1 (struct Out (make (do Len (link (1.0 . 2)))))
(native "libfftw3.so" "fftw_destroy_plan" NIL P)
(native "libfftw3.so" "fftw_free" NIL Out)
(native "libfftw3.so" "fftw_free" NIL In) ) ) )

Test:

(for R (fft '((1.0 0) (1.0 0) (1.0 0) (1.0 0) (0 0) (0 0) (0 0) (0 0)))
(tab (6 8)
(round (car R))
(round (cadr R)) ) )
Output:
 4.000   0.000
 1.000  -2.414
 0.000   0.000
 1.000  -0.414
 0.000   0.000
 1.000   0.414
 0.000   0.000
 1.000   2.414

[edit] PL/I

test: PROCEDURE OPTIONS (MAIN, REORDER); /* Derived from Fortran Rosetta Code */
 
/* In-place Cooley-Tukey FFT */
FFT: PROCEDURE (x) RECURSIVE;
DECLARE x(*) COMPLEX FLOAT (18);
DECLARE t COMPLEX FLOAT (18);
DECLARE ( N, Half_N ) FIXED BINARY (31);
DECLARE ( i, j ) FIXED BINARY (31);
DECLARE (even(*), odd(*)) CONTROLLED COMPLEX FLOAT (18);
DECLARE pi FLOAT (18) STATIC INITIAL ( 3.14159265358979323E0);
 
N = HBOUND(x);
 
if N <= 1 THEN return;
 
allocate odd((N+1)/2), even(N/2);
 
/* divide */
do j = 1 to N by 2; odd((j+1)/2) = x(j); end;
do j = 2 to N by 2; even(j/2) = x(j); end;
 
/* conquer */
call fft(odd);
call fft(even);
 
/* combine */
half_N = N/2;
do i=1 TO half_N;
t = exp(COMPLEX(0, -2*pi*(i-1)/N))*even(i);
x(i) = odd(i) + t;
x(i+half_N) = odd(i) - t;
end;
 
FREE odd, even;
 
END fft;
 
 
DECLARE data(8) COMPLEX FLOAT (18) STATIC INITIAL (
1, 1, 1, 1, 0, 0, 0, 0);
DECLARE ( i ) FIXED BINARY (31);
 
call fft(data);
 
do i=1 TO 8;
PUT SKIP LIST ( fixed(data(i), 25, 12) );
end;
 
END test;
Output:
    4.000000000000+0.000000000000I    
    1.000000000000-2.414213562373I    
    0.000000000000+0.000000000000I    
    1.000000000000-0.414213562373I    
    0.000000000000+0.000000000000I    
    0.999999999999+0.414213562373I    
    0.000000000000+0.000000000000I    
    0.999999999999+2.414213562373I

[edit] Prolog

Translation of: Python
Note: Similar algorithmically to the python example.
Works with: SWI Prolog version Version 6.2.6 by Jan Wielemaker, University of Amsterdam
:- dynamic twiddles/2.
%_______________________________________________________________
% Arithemetic for complex numbers; only the needed rules
add(cx(R1,I1),cx(R2,I2),cx(R,I)) :- R is R1+R2, I is I1+I2.
sub(cx(R1,I1),cx(R2,I2),cx(R,I)) :- R is R1-R2, I is I1-I2.
mul(cx(R1,I1),cx(R2,I2),cx(R,I)) :- R is R1*R2-I1*I2, I is R1*I2+R2*I1.
polar_cx(Mag, Theta, cx(R, I)) :- % Euler
R is Mag * cos(Theta), I is Mag * sin(Theta).
%___________________________________________________
% FFT Implementation. Note: K rdiv N is a rational number,
% making the lookup in dynamic database predicate twiddles/2 very
% efficient. Also, polar_cx/2 gets called only when necessary- in
% this case (N=8), exactly 3 times: (where Tf=1/4, 1/8, or 3/8).
tw(0,cx(1,0)) :- !. % Calculate e^(-2*pi*k/N)
tw(Tf, Cx) :- twiddles(Tf, Cx), !. % dynamic match?
tw(Tf, Cx) :- polar_cx(1.0, -2*pi*Tf, Cx), assert(twiddles(Tf, Cx)).
 
fftVals(N, Even, Odd, V0, V1) :- % solves all V0,V1 for N,Even,Odd
nth0(K,Even,E), nth0(K,Odd,O), Tf is K rdiv N, tw(Tf,Cx),
mul(Cx,O,M), add(E,M,V0), sub(E,M,V1).
 
split([],[],[]). % split [[a0,b0],[a1,b1],...] into [a0,a1,...] and [b0,b1,...]
split([[V0,V1]|T], [V0|T0], [V1|T1]) :- !, split(T, T0, T1).
 
fft([H], [H]).
fft([H|T], List) :-
length([H|T],N),
findall(Ve, (nth0(I,[H|T],Ve),I mod 2 =:= 0), EL), !, fft(EL, Even),
findall(Vo, (nth0(I,T,Vo),I mod 2 =:= 0),OL), !, fft(OL, Odd),
findall([V0,V1],fftVals(N,Even,Odd,V0,V1),FFTVals), % calc FFT
split(FFTVals,L0,L1), append(L0,L1,List).
%___________________________________________________
test :- D=[cx(1,0),cx(1,0),cx(1,0),cx(1,0),cx(0,0),cx(0,0),cx(0,0),cx(0,0)],
time(fft(D,DRes)), writef('fft=['), P is 10^3, !,
(member(cx(Ri,Ii), DRes), R is integer(Ri*P)/P, I is integer(Ii*P)/P,
write(R), (I>=0, write('+'),fail;write(I)), write('j, '),
fail; write(']'), nl).
 
Output:
 test.
% 681 inferences, 0.000 CPU in 0.001 seconds (0% CPU, Infinite Lips)
fft=[4+0j, 1-2.414j, 0+0j, 1-0.414j, 0+0j, 1+0.414j, 0+0j, 1+2.414j, ]
true.

[edit] Python

[edit] Python: Recursive

from cmath import exp, pi
 
def fft(x):
N = len(x)
if N <= 1: return x
even = fft2(x[0::2])
odd = fft2(x[1::2])
T= [exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]
return [even[k] + T[k] for k in xrange(N/2)] + \
[even[k] - T[k] for k in xrange(N/2)]
 
print( ' '.join("%5.3f" % abs(f)
for f in fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])) )
Output:
4.000 2.613 0.000 1.082 0.000 1.082 0.000 2.613

[edit] Python: Using module numpy

>>> from numpy.fft import fft
>>> from numpy import array
>>> a = array([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])
>>> print( ' '.join("%5.3f" % abs(f) for f in fft(a)) )
4.000 2.613 0.000 1.082 0.000 1.082 0.000 2.613

[edit] R

The function "fft" is readily available in R

fft(c(1,1,1,1,0,0,0,0))
Output:
4+0.000000i 1-2.414214i 0+0.000000i 1-0.414214i 0+0.000000i 1+0.414214i 0+0.000000i 1+2.414214i

[edit] Racket

 
#lang racket
(require math)
(array-fft (array #[1. 1. 1. 1. 0. 0. 0. 0.]))
 
Output:
(fcarray
 #[4.0+0.0i
   1.0-2.414213562373095i
   0.0+0.0i
   1.0-0.4142135623730949i
   0.0+0.0i
   0.9999999999999999+0.4142135623730949i
   0.0+0.0i
   0.9999999999999997+2.414213562373095i])

[edit] REXX

This REXX program is modeled after the Run BASIC version and is a radix-2 DIC (decimation-in-time) form of the
Cooley-Turkey FFT algorithm, and as such, this simplified form assumes that the number of data points is equal to an exact power of two.

Note that the REXX language doesn't have any higher math functions, so the functions COS and R2R
(normalization of radians) have been included here, as well as the constant pi .

This program also adds zero values if the number of data points in the list doesn't exactly equal to a power of two.
This is called zero-padding.

/*REXX pgm does a fast Fourier transform (FFT) on a set of complex nums.*/
numeric digits length( pi() ) - 1 /*limited by PI function result. */
arg data /*the ARG verb uppercases DATA */
if data='' then data=1 1 1 1 0 /*No data? Then use the default.*/
data=translate(data, 'J', "I") /*allow use of I as well as J */
size=words(data); pad=left('',6) /*PAD: for indenting/padding SAYs*/
do p=0 until 2**p>=size  ; end /* # args exactly a power of 2? */
do j=size+1 to 2**p;data=data 0; end /*add zeroes until a power of 2. */
size=words(data); ph=p%2; call hdr /*┌─────────────────────────────┐*/
/*│ Numbers in data can be in │*/
do j=0 for size /*│ 7 formats: real │*/
_=word(data,j+1) /*│ real,imag │*/
parse var _ #.1.j ',' #.2.j /*│ ,imag │*/
if right(#.1.j,1)=='J' then parse , /*│ nnnJ │*/
var #.1.j #2.j "J" @.1.j /*│ nnnj │*/
do m=1 for 2 /*omitted?*/ /*│ nnnI │*/
#.m.j=word(#.m.j 0, 1) /*│ nnni │*/
end /*m*/ /*└─────────────────────────────┘*/
 
say pad " FFT in " center(j+1,7) pad nice(#.1.j) nice(#.2.j,'i')
end /*j*/
 
say; say; tran=2*pi()/2**p;  !.=0
hp=2**p%2; counterA=2**(p-ph); pointer=counterA; doubler=1
 
do p-ph; halfpointer=pointer%2
 
do i=halfpointer by pointer to counterA-halfpointer
_=i-halfpointer;  !.i=!._+doubler
end /*i*/
 
doubler=doubler*2; pointer=halfpointer
end /*p-ph*/
 
do j=0 to 2**p%4; cmp.j=cos(j*tran); _m=hp-j; cmp._m=-cmp.j
_p=hp+j; cmp._p=-cmp.j
end /*j*/
 
counterB=2**ph
 
do i=0 for counterA; q=i *counterB
do j=0 for counterB; h=q+j; _=!.j*counterB+!.i; if _<=h then iterate
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._
end /*j*/ /* [↑] switch two sets of values*/
end /*i*/
 
double=1; do p  ; w=hp%double
do k=0 for double; lb=w*k  ; lh=lb+2**p%4
do j=0 for w  ; a=j*double*2+k  ; b=a+double
r=#.1.a; i=#.2.a ; c1=cmp.lb*#.1.b ; c4=cmp.lb*#.2.b
c2=cmp.lh*#.2.b ; c3=cmp.lh*#.1.b
#.1.a=r+c1-c2  ; #.2.a=i+c3+c4
#.1.b=r-c1+c2  ; #.2.b=i-c3-c4
end /*j*/
end /*k*/
double=double+double
end /*p*/
call hdr
do i=0 for size
say pad " FFT out " center(i+1,7) pad nice(#.1.i) nice(#.2.i,'j')
end /*i*/
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────HDR subroutine──────────────────────*/
hdr: _='───data─── num real-part imaginary-part'; say pad _
say pad translate(_, " "copies('═',256), " "xrange()); return
/*──────────────────────────────────PI subroutine───────────────────────────────────*/
pi: return , /*add more digs if NUMERIC DIGITS > 85. */
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628
/*──────────────────────────────────R2R subroutine──────────────────────*/
r2r: return arg(1) // (2*pi()) /*reduce radians to unit circle. */
/*──────────────────────────────────COS subroutine──────────────────────*/
cos: procedure; parse arg x; x=r2r(x); return .sincos(1,1,-1)
.sincos: parse arg z,_,i; x=x*x; p=z
do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end
return z
/*──────────────────────────────────NICE subroutine─────────────────────*/
nice: procedure; parse arg x,j /*makes complex nums look nicer. */
numeric digits digits()%10; nz='1e-'digits() /*show ≈10% of DIGITS.*/
if abs(x)<nz then x=0; x=x/1; if x=0 & j\=='' then return ''
x=format(x,,digits()); if pos('.',x)\==0 then x=strip(x,'T',0)
x=strip(x,,'.'); if x>=0 then x=' '||x; return left(x||j,digits()+4)
Output:
when using the default input of:   1 1 1 1 0
       ───data───   num       real-part      imaginary-part
       ══════════   ═══       ═════════      ══════════════
        FFT   in     1            1
        FFT   in     2            1
        FFT   in     3            1
        FFT   in     4            1
        FFT   in     5            0
        FFT   in     6            0
        FFT   in     7            0
        FFT   in     8            0


       ───data───   num       real-part      imaginary-part
       ══════════   ═══       ═════════      ══════════════
        FFT  out     1            4
        FFT  out     2            1           -2.4142136j
        FFT  out     3            0
        FFT  out     4            1           -0.41421356j
        FFT  out     5            0
        FFT  out     6            1            0.41421356j
        FFT  out     7            0
        FFT  out     8            1            2.4142136j

[edit] Ruby

def fft(vec)
return vec if vec.size <= 1
evens_odds = vec.partition.with_index{|_,i| i.even?}
evens, odds = evens_odds.map{|even_odd| fft(even_odd)*2}
evens.zip(odds).map.with_index do |(even, odd),i|
even + odd * Math::E ** Complex(0, 2 * Math::PI * (-i)/ vec.size)
end
end
 
fft([1,1,1,1,0,0,0,0]).each{|c| p c}
 
Output:
(4+0i)
(1.0-2.414213562373095i)
(-1.2246467991473532e-16-1.2246467991473532e-16i)
(1.0-0.4142135623730949i)
(0.0-2.4492935982947064e-16i)
(0.9999999999999998+0.41421356237309515i)
(1.2246467991473532e-16-1.224646799147353e-16i)
(0.9999999999999993+2.414213562373095i)


[edit] Run BASIC

cnt  = 8
sig = int(log(cnt) /log(2) +0.9999)
 
pi = 3.14159265
real1 = 2^sig
 
real = real1 -1
real2 = int(real1 / 2)
real4 = int(real1 / 4)
real3 = real4 +real2
 
dim rel(real1)
dim img(real1)
dim cmp(real3)
 
for i = 0 to cnt -1
read rel(i)
read img(i)
next i
 
data 1,0, 1,0, 1,0, 1,0, 0,0, 0,0, 0,0, 0,0
 
sig2 = int(sig / 2)
sig1 = sig -sig2
cnt1 = 2^sig1
cnt2 = 2^sig2
 
dim v(cnt1 -1)
v(0) = 0
dv = 1
ptr = cnt1
 
for j = 1 to sig1
hlfPtr = int(ptr / 2)
pt = cnt1 - hlfPtr
for i = hlfPtr to pt step ptr
v(i) = v(i -hlfPtr) + dv
next i
dv = dv + dv
ptr = hlfPtr
next j
 
k = 2 *pi /real1
 
for x = 0 to real4
cmp(x) = cos(k *x)
cmp(real2 - x) = 0 - cmp(x)
cmp(real2 + x) = 0 - cmp(x)
next x
 
print "fft: bit reversal"
 
for i = 0 to cnt1 -1
ip = i *cnt2
for j = 0 to cnt2 -1
h = ip +j
g = v(j) *cnt2 +v(i)
if g >h then
temp = rel(g)
rel(g) = rel(h)
rel(h) = temp
temp = img(g)
img(g) = img(h)
img(h) = temp
end if
next j
next i
 
t = 1
for stage = 1 to sig
print " stage:- "; stage
d = int(real2 / t)
for ii = 0 to t -1
l = d *ii
ls = l +real4
for i = 0 to d -1
a = 2 *i *t +ii
b = a +t
f1 = rel(a)
f2 = img(a)
cnt1 = cmp(l) *rel(b)
cnt2 = cmp(ls) *img(b)
cnt3 = cmp(ls) *rel(b)
cnt4 = cmp(l) *img(b)
rel(a) = f1 + cnt1 - cnt2
img(a) = f2 + cnt3 + cnt4
rel(b) = f1 - cnt1 + cnt2
img(b) = f2 - cnt3 - cnt4
next i
next ii
t = t +t
next stage
 
print " Num real imag"
for i = 0 to real
if abs(rel(i)) <10^-5 then rel(i) = 0
if abs(img(i)) <10^-5 then img(i) = 0
print " "; i;" ";using("##.#",rel(i));" ";img(i)
next i
end
  Num   real   imag
   0    4.0    0
   1    1.0    -2.41421356
   2    0.0    0
   3    1.0    -0.414213565
   4    0.0    0
   5    1.0    0.414213562
   6    0.0    0
   7    1.0    2.41421356

[edit] Scala

Library: Scala
Works with: Scala version 2.10.4

Imports and Complex arithmetic:

import scala.math.{ Pi, cos, sin, cosh, sinh, abs }
 
case class Complex(re: Double, im: Double) {
def +(x: Complex): Complex = Complex(re + x.re, im + x.im)
def -(x: Complex): Complex = Complex(re - x.re, im - x.im)
def *(x: Double): Complex = Complex(re * x, im * x)
def *(x: Complex): Complex = Complex(re * x.re - im * x.im, re * x.im + im * x.re)
def /(x: Double): Complex = Complex(re / x, im / x)
 
override def toString(): String = {
val a = "%1.3f" format re
val b = "%1.3f" format abs(im)
(a,b) match {
case (_, "0.000") => a
case ("0.000", _) => b + "i"
case (_, _) if im > 0 => a + " + " + b + "i"
case (_, _) => a + " - " + b + "i"
}
}
}
 
def exp(c: Complex) : Complex = {
val r = (cosh(c.re) + sinh(c.re))
Complex(cos(c.im), sin(c.im)) * r
}

The FFT definition itself:

def _fft(cSeq: Seq[Complex], direction: Complex, scalar: Int): Seq[Complex] = {
if (cSeq.length == 1) {
return cSeq
}
val n = cSeq.length
assume(n % 2 == 0, "The Cooley-Tukey FFT algorithm only works when the length of the input is even.")
 
val evenOddPairs = cSeq.grouped(2).toSeq
val evens = _fft(evenOddPairs map (_(0)), direction, scalar)
val odds = _fft(evenOddPairs map (_(1)), direction, scalar)
 
def leftRightPair(k: Int): Pair[Complex, Complex] = {
val base = evens(k) / scalar
val offset = exp(direction * (Pi * k / n)) * odds(k) / scalar
(base + offset, base - offset)
}
 
val pairs = (0 until n/2) map leftRightPair
val left = pairs map (_._1)
val right = pairs map (_._2)
left ++ right
}
 
def fft(cSeq: Seq[Complex]): Seq[Complex] = _fft(cSeq, Complex(0, 2), 1)
def rfft(cSeq: Seq[Complex]): Seq[Complex] = _fft(cSeq, Complex(0, -2), 2)

Usage:

val data = Seq(Complex(1,0), Complex(1,0), Complex(1,0), Complex(1,0), 
Complex(0,0), Complex(0,2), Complex(0,0), Complex(0,0))
 
println(fft(data))
println(rfft(fft(data)))
Output:
Vector(4.000 + 2.000i, 2.414 + 1.000i, -2.000, 2.414 + 1.828i, 2.000i, -0.414 + 1.000i, 2.000, -0.414 - 3.828i)
Vector(1.000, 1.000, 1.000, 1.000, 0.000, 2.000i, 0.000, 0.000)

[edit] Scilab

Scilab has a builtin FFT function.

fft([1,1,1,1,0,0,0,0]')

[edit] Tcl

Library: Tcllib (Package: math::constants)
Library: Tcllib (Package: math::fourier)
package require math::constants
package require math::fourier
 
math::constants::constants pi
# Helper functions
proc wave {samples cycles} {
global pi
set wave {}
set factor [expr {2*$pi * $cycles / $samples}]
for {set i 0} {$i < $samples} {incr i} {
lappend wave [expr {sin($factor * $i)}]
}
return $wave
}
proc printwave {waveName {format "%7.3f"}} {
upvar 1 $waveName wave
set out [format "%-6s" ${waveName}:]
foreach value $wave {
append out [format $format $value]
}
puts $out
}
proc waveMagnitude {wave} {
set out {}
foreach value $wave {
lassign $value re im
lappend out [expr {hypot($re, $im)}]
}
return $out
}
 
set wave [wave 16 3]
printwave wave
# Uses FFT if input length is power of 2, and a less efficient algorithm otherwise
set fft [math::fourier::dft $wave]
# Convert to magnitudes for printing
set fft2 [waveMagnitude $fft]
printwave fft2
Output:
wave:   0.000  0.924  0.707 -0.383 -1.000 -0.383  0.707  0.924  0.000 -0.924 -0.707  0.383  1.000  0.383 -0.707 -0.924
fft2:   0.000  0.000  0.000  8.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  8.000  0.000  0.000

[edit] Ursala

The fftw library is callable from Ursala using the syntax ..u_fw_dft for a one dimensional forward discrete Fourier transform operating on a list of complex numbers. Ordinarily the results are scaled so that the forward and reverse transforms are inverses of each other, but additional scaling can be performed as shown below to conform to convention.

#import nat
#import flo
 
f = <1+0j,1+0j,1+0j,1+0j,0+0j,0+0j,0+0j,0+0j> # complex sequence of 4 1's and 4 0's
 
g = c..mul^*D(sqrt+ float+ length,..u_fw_dft) f # its fft
 
#cast %jLW
 
t = (f,g)
Output:
(
   <
      1.000e+00+0.000e+00j,
      1.000e+00+0.000e+00j,
      1.000e+00+0.000e+00j,
      1.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j,
      0.000e+00+0.000e+00j>,
   <
      4.000e+00+0.000e+00j,
      1.000e+00-2.414e+00j,
      0.000e+00+0.000e+00j,
      1.000e+00-4.142e-01j,
      0.000e+00+0.000e+00j,
      1.000e+00+4.142e-01j,
      0.000e+00+0.000e+00j,
      1.000e+00+2.414e+00j>)
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox