Fast Fourier transform: Difference between revisions

Content added Content deleted

Inline

Revision as of 20:20, 1 November 2015

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.

Ada

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

<lang Ada> 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;

</lang>

<lang Ada> 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; </lang>

Example:

<lang Ada> 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; </lang>

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)

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<lang algol68>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

);</lang>File: test.Fast_Fourier_transform.a68<lang algol68>#!/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$
 ))

)</lang>

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,

BBC BASIC

Works with: BBC BASIC for Windows

<lang bbcbasic> @% = &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

</lang>

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

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.<lang C>

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; }

</lang>

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)

C++

<lang cpp>#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, divide-and-conquer) // Higher memory requirements and redundancy although more intuitive 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;
   }

}

// Cooley-Tukey FFT (in-place, breadth-first, decimation-in-frequency) // Better optimized but less intuitive void fft(CArray &x) { // DFT unsigned int N = x.size(), k = N, n; double thetaT = 3.14159265358979323846264338328L / N; Complex phiT = Complex(cos(thetaT), sin(thetaT)), T; while (k > 1) { n = k; k >>= 1; phiT = phiT * phiT; T = 1.0L; for (unsigned int l = 0; l < k; l++) { for (unsigned int a = l; a < N; a += n) { unsigned int b = a + k; Complex t = x[a] - x[b]; x[a] += x[b]; x[b] = t * T; } T *= phiT; } } // Decimate unsigned int m = (unsigned int)log2(N); for (unsigned int a = 0; a < N; a++) { unsigned int b = a; // Reverse bits b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1)); b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2)); b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4)); b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8)); b = ((b >> 16) | (b << 16)) >> (32 - m); if (b > a) { Complex t = x[a]; x[a] = x[b]; x[b] = t; } } //// Normalize (This section make it not working correctly) //Complex f = 1.0 / sqrt(N); //for (unsigned int i = 0; i < N; i++) // x[i] *= f; }

// 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;

}</lang>

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)

Common Lisp

Translation of: Python

<lang lisp> (defun fft (x)

 (if (<= (length x) 1) x 
   (let*
      (
        (even (fft (loop for i from 0 below (length x) by 2 collect (nth i x))))
        (odd (fft (loop for i from 1 below (length x) by 2 collect (nth i x))))
        (aux (loop for k from 0 below (/ (length x) 2) collect (* (exp (/ (* (complex 0 -2) pi k ) (length x))) (nth k odd))))
      )
      (append (mapcar #'+ even aux) (mapcar #'- even aux))
      )
   )

)

(mapcar (lambda (x) (format t "~a~&" x)) (fft '(1 1 1 1 0 0 0 0)))

</lang>

Output:

#C(4.0d0 0.0d0)
#C(1.0d0 -2.414213562373095d0)
#C(0.0d0 0.0d0)
#C(1.0d0 -0.4142135623730949d0)
#C(0.0d0 0.0d0)
#C(0.9999999999999999d0 0.4142135623730949d0)
#C(0.0d0 0.0d0)
#C(0.9999999999999997d0 2.414213562373095d0)

D

Standard Version

<lang d>void main() {

   import std.stdio, std.numeric;

   [1.0, 1, 1, 1, 0, 0, 0, 0].fft.writeln;

}</lang>

Output:

[4+0i, 1-2.41421i, 0-0i, 1-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

creals Version

Built-in complex numbers will be deprecated. <lang d>import std.stdio, std.algorithm, std.range, std.math;

const(creal)[] fft(in creal[] x) pure /*nothrow*/ @safe {

   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() @safe {

   [1.0L+0i, 1, 1, 1, 0, 0, 0, 0].fft.writeln;

}</lang>

Output:

[4+0i, 1+-2.41421i, 0+0i, 1+-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

Phobos Complex Version

<lang d>import std.stdio, std.algorithm, std.range, std.math, std.complex;

auto fft(T)(in T[] x) pure /*nothrow @safe*/ {

   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;

}</lang>

Output:

[4+0i, 1-2.41421i, 0+0i, 1-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]

ERRE

<lang ERRE> PROGRAM FFT

CONST CNT=8

!$DYNAMIC DIM REL[0],IMG[0],CMP[0],V[0]

BEGIN

  SIG=INT(LOG(CNT)/LOG(2)+0.9999)
  REAL1=2^SIG

  REAL=REAL1-1
  REAL2=INT(REAL1/2)
  REAL4=INT(REAL1/4)
  REAL3=REAL4+REAL2

!$DIM REL[REAL1],IMG[REAL1],CMP[REAL3]

FOR I=0 TO CNT-1 DO

  READ(REL[I],IMG[I])

END FOR

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 DO

 HLFPTR=INT(PTR/2)
 PT=CNT1-HLFPTR
 FOR I=HLFPTR TO PT STEP PTR DO
   V[I]=V[I-HLFPTR]+DV
 END FOR
 DV=2*DV
 PTR=HLFPTR

END FOR

K=2*π/REAL1

FOR X=0 TO REAL4 DO

  CMP[X]=COS(K*X)
  CMP[REAL2-X]=-CMP[X]
  CMP[REAL2+X]=-CMP[X]

END FOR

PRINT("FFT: BIT REVERSAL")

FOR I=0 TO CNT1-1 DO

 IP=I*CNT2
 FOR J=0 TO CNT2-1 DO
   H=IP+J
   G=V[J]*CNT2+V[I]
   IF G>H THEN
      SWAP(REL[G],REL[H])
      SWAP(IMG[G],IMG[H])
   END IF
 END FOR

END FOR

T=1 FOR STAGE=1 TO SIG DO

 PRINT("STAGE:";STAGE)
 D=INT(REAL2/T)
 FOR II=0 TO T-1 DO
    L=D*II
    LS=L+REAL4
    FOR I=0 TO D-1 DO
      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
    END FOR
 END FOR
 T=2*T

END FOR

PRINT("NUM REAL IMAG") FOR I=0 TO REAL DO

   IF ABS(REL[I])<1E-5 THEN REL[I]=0 END IF
   IF ABS(IMG[I])<1E-5 THEN IMG[I]=0 END IF
   PRINT(I;"";)
   WRITE("##.###### ##.######";REL[I];IMG[I])

END FOR END PROGRAM </lang>

Output:

FFT: BIT REVERSAL
STAGE: 1
STAGE: 2
STAGE: 3
NUM REAL     IMAG
 0  4.000000  0.000000
 1  1.000000 -2.414214
 2  0.000000  0.000000
 3  1.000000 -0.414214
 4  0.000000  0.000000
 5  1.000000  0.414214
 6  0.000000  0.000000
 7  1.000000  2.414214

Factor

<lang 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 }

} </lang>

Fortran

<lang 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</lang>

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 )

GAP

<lang 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 ]</lang>

Go

<lang 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)
   }

}</lang>

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)

Haskell

<lang 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]</lang>

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

J

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

<lang j>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@#</lang>

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):

Here is a representation of an example which appears in some of the other implementations, here: <lang J> 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</lang>

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

JavaScript

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

<lang javascript>/* 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])) );</lang> Very very basic Complex number that provides only the components required by the code above. <lang javascript>/* 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'); }</lang>

jq

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

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) ];</lang>

FFT

<lang jq>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;</lang>

Example: <lang jq>[1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0] | fft </lang>

Output:

[[4,-0],[1,-2.414213562373095],
 [0,0],[1,-0.4142135623730949],
 [0,0],[0.9999999999999999,0.4142135623730949],
 [0,0],[0.9999999999999997,2.414213562373095]]

Julia

Julia has a built-in FFT function: <lang julia>fft([1,1,1,1,0,0,0,0])</lang>

Output:

<lang julia>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</lang>

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

</lang>

   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

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.

<lang Maple> with( DiscreteTransforms ):

FourierTransform( <1,1,1,1,0,0,0,0>, normalization=none ); </lang>

                         [       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. <lang Maple> v := Vector( [1,1,1,1,0,0,0,0], datatype=complex[8] ):

FourierTransform( v, normalization=none, inplace ):

v; </lang>

                         [       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 ]

<lang Maple> InverseFourierTransform( v, normalization=full, inplace ):

v; </lang>

                       [          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          ]

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.

<lang Mathematica> Fourier[{1,1,1,1,0,0,0,0}, FourierParameters->{1,-1}] </lang>

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}

MATLAB / Octave

Matlab/Octave have a builtin FFT function.

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

Maxima

Nim

Translation of: Python

<lang nim>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)</lang>

Output:

4.000
2.613
-0.000
1.082
-0.000
1.082
-0.000
2.613

OCaml

This is a simple implementation of the Cooley-Tukey pseudo-code <lang OCaml>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)</lang>

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)

PARI/GP

Naive implementation, using the same testcase as Ada: <lang parigp>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])</lang>

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]

Perl

Translation of: Perl 6

<lang Perl>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);</lang>

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

Perl 6

Works with: rakudo version 2014.12-39-g7ea485c built on MoarVM version 2014.12

<lang perl6>sub fft {

   return @_ if @_ == 1;
   my @evn = fft( @_[0, 2 ... *] );
   my @odd = fft( @_[1, 3 ... *] ) Z*
   map &cis, (0, 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");</lang>

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

PicoLisp

Works with: PicoLisp version 3.1.0.3

<lang PicoLisp># 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) ) ) )</lang>

Test: <lang PicoLisp>(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)) ) )</lang>

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

PL/I

<lang 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;</lang>

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

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

<lang prolog>:- 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). </lang>

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.

Python

Python: Recursive

<lang python>from cmath import exp, pi

def fft(x):

   N = len(x)
   if N <= 1: return x
   even = fft(x[0::2])
   odd =  fft(x[1::2])
   T= [exp(-2j*pi*k/N)*odd[k] for k in range(N//2)]
   return [even[k] + T[k] for k in range(N//2)] + \
          [even[k] - T[k] for k in range(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])) )</lang>

Output:

4.000 2.613 0.000 1.082 0.000 1.082 0.000 2.613

Python: Using module numpy

<lang python>>>> 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</lang>

R

The function "fft" is readily available in R <lang R>fft(c(1,1,1,1,0,0,0,0))</lang>

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

Racket

lang racket

(require math) (array-fft (array #[1. 1. 1. 1. 0. 0. 0. 0.])) </lang>

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])

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, such as the functions COS and R2R
(cosine and reduce radians to a unit circle).

A normalization of radians function has been included here, as well as the constant pi.

This REXX 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 known as zero-padding. <lang rexx>/*REXX pgm performs a fast Fourier transform (FFT) on a set of complex numbers*/ numeric digits length( pi() ) - 1 /*limited by the PI function result. */ arg data /*ARG verb uppercases the DATA from CL.*/ if data= then data=1 1 1 1 0 /*Not specified? Then use the default.*/ size=words(data); pad=left(,6) /*PAD: for indenting and padding SAYs.*/

 do p=0  until  2**p>=size      ; end /*number of args exactly a power of 2? */
 do j=size+1 to 2**p;data=data 0; end /*add zeroes to DATA 'til a power of 2.*/

size=words(data); ph=p%2; call hdr /*┌─────────────────────────────┐*/

            /* [↓] TRANSLATE allows I&J*/   /*│ Numbers in  data  can be in │*/
    do j=0  for size                        /*│ seven formats:    real      │*/
    _=translate(word(data,j+1), 'J', "I")   /*│                   real,imag │*/
    parse  var  _   #.1.j  $ 1 ',' #.2.j  /*│                       ,imag │*/
    if $=='J'  then parse var #.1.j #2.j ,  /*│                        nnnJ │*/
                              "J"  #.1.j    /*│                        nnnj │*/
      do m=1  for  2; #.m.j=word(#.m.j 0,1) /*│                        nnnI │*/
      end   /*m*/   /* [↑] ommited part?*/  /*│                        nnni │*/
                                            /*└─────────────────────────────┘*/
    say pad " FFT   in "    center(j+1,7)    pad   fmt(#.1.j)    fmt(#.2.j,'i')
    end     /*j*/

say tran=pi()*2/2**p; !.=0; hp=2**p%2; A=2**(p-ph); ptr=A; dbl=1 say

   do p-ph;        halfPtr=ptr%2
              do i=halfPtr  by ptr  to A-halfPtr;    _=i-halfPtr;   !.i=!._+dbl
              end   /*i*/
   dbl=dbl*2;               ptr=halfPtr
   end   /*p-ph*/

            do j=0  to 2**p%4;  cmp.j=cos(j*tran);    _=hp - j;   cmp._= -cmp.j
                                                      _=hp + j;   cmp._= -cmp.j
            end  /*j*/

B=2**ph

 do i=0     for A;            q=i * B
     do j=0 for B;   h=q+j;   _=!.j*B+!.i;  if _<=h  then iterate
     parse value  #.1._  #.1.h  #.2._  #.2.h   with  #.1.h  #.1._  #.2.h  #.2._
     end   /*j*/                             /* [↑]  swap two sets of values.*/
 end       /*i*/

dbl=1; do p ; w=hp % dbl

              do k=0   for dbl   ;   Lb=w * k            ;    Lh=Lb + 2**p % 4
                do j=0 for w     ;    a=j * dbl * 2 + k  ;     b= a + dbl
                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*/
            dbl=dbl+dbl
            end         /*p*/

call hdr

        do i=0  for size
        say pad  " FFT  out "   center(i+1,7)  pad  fmt(#.1.i)   fmt(#.2.i,'j')
        end   /*i*/

exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ cos: procedure; parse arg x; q=r2r(x)**2; z=1; _=1; p=1

    do k=2 by 2; _=-_*q/(k*(k-1)); z=z+_; if z=p then leave; p=z; end; return z

/*────────────────────────────────────────────────────────────────────────────*/ hdr: _='───data─── num real─part imaginary─part'; say pad _

    say pad  translate(_, " "copies('═',256), " "xrange());              return

/*────────────────────────────────────────────────────────────────────────────*/ pi: return 3.141592653589793238462643383279502884197169399375105820974944592308 /*────────────────────────────────────────────────────────────────────────────*/ fmt: procedure; parse arg y,j /*transforms complex numbers look fmtr*/

    numeric digits digits()%8;        nz='1e-'digits()  /*show ≈10% of DIGITS*/
    if abs(y)<nz  then y=0;           y=y/1;    if y=0 & j\==  then return 
    y=format(y,,digits());            if pos(.,y)\==0  then y=strip(y,'T',0)
    y=strip(y,,.);   if y>=0  then y=' 'y;        return left(y||j,digits()+4)

/*────────────────────────────────────────────────────────────────────────────*/ r2r: return arg(1) // (pi()*2) /*reduce the radians to a unit circle. */</lang> 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

Ruby

<lang 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| puts "%9.6f %+9.6fi" % c.rect}</lang>

Output:

 4.000000 +0.000000i
 1.000000 -2.414214i
-0.000000 -0.000000i
 1.000000 -0.414214i
 0.000000 -0.000000i
 1.000000 +0.414214i
 0.000000 -0.000000i
 1.000000 +2.414214i

Run BASIC

<lang runbasic>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</lang>

  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

Scala

Library: Scala

Works with: Scala version 2.10.4

Imports and Complex arithmetic: <lang Scala>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

}</lang>

The FFT definition itself: <lang Scala>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)</lang>

Usage: <lang Scala>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)))</lang>

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)

Scilab

Scilab has a builtin FFT function.

Sidef

Translation of: Perl

<lang ruby>func fft(arr) {

   arr.len == 1 && return arr;

   var evn = fft([arr.@[arr.range.grep { .is_even }]]);
   var odd = fft([arr.@[arr.range.grep { .is_odd  }]]);
   var twd = (Complex(0, 2) * Math::PI / arr.len);

   odd.range.map {|n| odd[n] *= exp(twd * n)};
   (evn »+« odd) + (evn »-« odd);

}

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

Tcl

Library: Tcllib (Package: math::constants)

Library: Tcllib (Package: math::fourier)

<lang tcl>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</lang>

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

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. <lang ursala>#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)</lang>

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