Fast Fourier transform: Difference between revisions

From Rosetta Code
Content added Content deleted
Line 114: Line 114:
INT lwb = (lwb OF slice|(INT lwb): lwb|LWB array);
INT lwb = (lwb OF slice|(INT lwb): lwb|LWB array);
INT upb = (upb OF slice|(INT upb): upb|UPB array);
INT upb = (upb OF slice|(INT upb): upb|UPB array);
(FALSE|printf(($"ELEM array in: "$, slice fmt, lwb, by, upb, $f(compl fmt)x$, array,$" => "$)));


[(upb-lwb)OVERby+1]VALUE out;
[(upb-lwb)%by+1]VALUE out;
INT upb out := 0;
INT upb out := 0;
FOR index FROM lwb BY by TO upb DO
FOR index FROM lwb BY by TO upb DO
out[upb out+:=1] := array[index] OD;
out[upb out+:=1] := array[index] OD;
(FALSE|printf(($"out: "$, slice fmt, lwb, by, upb, $f(compl fmt)x$, out,$l$)));
out
out
);
);
Line 128: Line 126:
ELSE
ELSE
INT n = UPB in - LWB in + 1;
INT n = UPB in - LWB in + 1;

[LWB in: UPB in]VALUE out;
[LWB in: UPB in]VALUE out;


[]VALUE even = fft(SLICE(LWB in, 2,EMPTY) ELEM in),
[]VALUE even = fft(SLICE(EMPTY, 2,EMPTY) ELEM in),
odd = fft(SLICE(LWB in+1,2,EMPTY) ELEM in);
odd = fft(SLICE(LWB in+1,2,EMPTY) ELEM in);
INT mid in = ( UPB in + LWB in ) OVER 2;
INT mid in = (UPB in + LWB in) % 2;
FOR i FROM LWB in TO mid in DO
FOR i FROM LWB in TO mid in DO
INT k = i - LWB in;
INT k = i - LWB in;
Line 143: Line 141:
FI
FI
);
);

FORMAT real fmt = $g(0,real width OVER 3-2)$;
FORMAT real fmt = $g(0,real width % 3-2)$;
FORMAT compl fmt = $f(real fmt)"⊥"f(real fmt)$;
FORMAT compl fmt = $f(real fmt)"⊥"f(real fmt)$;
FORMAT array compl fmt = $f(compl fmt)", "$;
FORMAT array compl fmt = $f(compl fmt)", "$;
Line 150: Line 148:
printf((
printf((
$f(array compl fmt)$, fft((1, 1, 1, 1, 0, 0, 0, 0)), $l$,
$f(array compl fmt)$, fft((1, 1, 1, 1, 0, 0, 0, 0)), $l$,
$f(array compl fmt)$,
$f(array compl fmt)$,
fft((0, 0.924, 0.707,-0.383,-1,-0.383, 0.707, 0.924,
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)), $l$
0,-0.924,-0.707, 0.383, 1, 0.383,-0.707,-0.924)), $l$
))</lang>'''Output:'''
))</lang>'''Output:'''

Revision as of 08:07, 4 December 2011

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.

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 - no extension to language.
Works with: ALGOL 68G version Any - tested with release algol68g-2.3.5.

File: Fast_Fourier_transform.a68<lang algol68>#!/usr/local/bin/a68g --script #

  1. -*- coding: utf-8 -*- #

MODE VALUE = COMPL;

MODE OPTINT = UNION(INT, VOID), 

    SLICE = STRUCT(OPTINT lwb, by, upb);

FORMAT slice fmt = $"["g(0)":"g(0)":"g(0)"] "$;

OP ELEM = (SLICE slice, []VALUE array)[]VALUE: ( 

   INT by = (by OF slice|(INT by): by|1);
   INT lwb = (lwb OF slice|(INT lwb): lwb|LWB array);
   INT upb = (upb OF slice|(INT upb): upb|UPB array);

   [(upb-lwb)%by+1]VALUE out;
   INT upb out := 0;
   FOR index FROM lwb BY by TO upb DO
     out[upb out+:=1] := array[index] OD;
   out

);

PROC fft = ([]VALUE in)[]VALUE: ( 

   IF LWB in >= UPB in THEN in
   ELSE
       INT n = UPB in - LWB in + 1;

       [LWB in: UPB in]VALUE out;

       []VALUE even = fft(SLICE(EMPTY,   2,EMPTY) ELEM in),
               odd =  fft(SLICE(LWB in+1,2,EMPTY) ELEM in);
       INT mid in = (UPB in + LWB in) % 2;
       FOR i FROM LWB in TO mid in DO
           INT k = i - LWB in;
           COMPL exp = complex exp(-(0I2)*pi*k/n)*odd[i];
           out[i]          := even[i] + exp;
           out[mid in + i] := even[i] - exp
       OD;
       out
   FI

);

FORMAT real fmt = $g(0,real width % 3-2)$; FORMAT compl fmt = $f(real fmt)"⊥"f(real fmt)$; FORMAT array compl fmt = $f(compl fmt)", "$;

printf(( 

 $f(array compl fmt)$, fft((1, 1, 1, 1, 0, 0, 0, 0)), $l$,
 $f(array compl fmt)$,
   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)), $l$

))</lang>Output: 

4.000⊥.000, 1.000⊥-2.414, .000⊥.000, 1.000⊥-.414, .000⊥.000, 1.000⊥.414, .000⊥.000, 1.000⊥2.414, 
.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,

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>

  1. include <math.h>
  2. 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); }

int main() { PI = atan2(1, 1) * 4; cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};

void show(char * s) { printf(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])); }

show("Data: "); fft(buf, 8); show("\nFFT : ");

return 0; }</lang>Output<lang>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)</lang>

C++

<lang cpp>#include <complex>

  1. include <iostream>
  2. 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;
   }

}

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

   fft(data);

   for (int i = 0; i < 8; ++i)
   {
       std::cout << data[i] << "\n";
   }
   return 0;

}</lang> Output: 

(4,0)
(1,-2.41421)
(0,0)
(1,-0.414214)
(0,0)
(1,0.414214)
(0,0)
(1,2.41421)

D

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

auto fft(creal[] x) { 

 int N = x.length;
 if (N <= 1) return x;
 auto ev = fft(array(stride(x, 2)));
 auto od = fft(array(stride(x[1 .. $], 2)));
 auto l = map!((k){return ev[k]+expi(-2*PI*k/N)*od[k];})(iota(N/2));
 auto r = map!((k){return ev[k]-expi(-2*PI*k/N)*od[k];})(iota(N/2));
 return array(chain(l, r));

}

void main() { 

 writeln(fft([1.0L+0i, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]));

}</lang>

Output:

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

GAP

<lang gap># Here an implementation with no optimization (O(n^2)).

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

  1. [ 4, 1-E(8)-E(8)^2-E(8)^3, 0, 1-E(8)+E(8)^2-E(8)^3,
  2. 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, 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 Control.Monad import Data.Array import Data.Complex

-- Cooley-Tukey fft [] = [] fft [x] = [x] fft xs = zipWith (+) ys ts ++ zipWith (-) ys ts 

   where n = length xs
         xs' = listArray (1,n) xs
         ys = fft $ map (xs'!) [1,3..n]
         zs = fft $ map (xs'!) [2,4..n]
         ts = zipWith (\z k -> (exp' k n) * z) zs [0..(n `div` 2)-1]
         exp' k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)

main = 

   do let res = fft [1,1,1,1,0,0,0,0]
      forM_ res print</lang>

And the 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:

<lang j> require'printf' 

  fmt =: [:, sprintf~&'%7.3f'"0
  ('wave:',:'fft:'),.fmt"1 (,: |@fft) 1 o. 2p1*3r16 * i.16

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

Note that sprintf does not support complex arguments, so we only display the magnitude of the fft here. 

  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

Liberty BASIC

<lang lb> 

   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


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}

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 6

<lang perl6>sub fft { 

   return @_ if @_ == 1;
   my @evn = fft( @_[0,2...^* >= @_] );
   my @odd = fft( @_[1,3...^* >= @_] );
   my $twd = 2i * pi / @_; # twiddle factor
   @odd  »*=« (^@odd).map(* * $twd)».exp;
   return @evn »+« @odd, @evn »-« @odd;

}

my @seq = ^16; my $cycles = 3; my @wave = (@seq »*» (2*pi / @seq * $cycles))».sin; 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

Python

<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])
   return [even[k] + exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)] + \
          [even[k] - exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]

print fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])</lang>

Output: 

[(4+0j), (1-2.4142135623730949j), 0j, (1-0.41421356237309492j), 0j, (0.99999999999999989+0.41421356237309492j), 0j, (0.99999999999999967+2.4142135623730949j)]

Using module numpy

<lang python>>>> from numpy.fft import fft >>> from numpy import array >>> a = array((0.0, 0.924, 0.707, -0.383, -1.0, -0.383, 0.707, 0.924, 0.0, -0.924, -0.707, 0.383, 1.0, 0.383, -0.707, -0.924)) >>> print( ' '.join("%5.3f" % abs(f) for f in fft(a)) ) 0.000 0.001 0.000 8.001 0.000 0.001 0.000 0.001 0.000 0.001 0.000 0.001 0.000 8.001 0.000 0.001</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

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

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

  1. Uses FFT if input length is power of 2, and a less efficient algorithm otherwise

set fft [math::fourier::dft $wave]

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

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

  1. 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>)
Retrieved from "https://rosettacode.org/wiki/Fast_Fourier_transform?oldid=171453"