Fast Fourier transform: Difference between revisions

From Rosetta Code
Content added Content deleted
(Ada solution added)
Line 1: Line 1:
{{task}}
{{task}}
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. [http://en.wikipedia.org/wiki/Cooley–Tukey_FFT_algorithm Wikipedia] has pseudocode for that. Further optimizations are possible but not required.
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. [http://en.wikipedia.org/wiki/Cooley–Tukey_FFT_algorithm Wikipedia] has pseudocode for that. Further optimizations are possible but not required.
=={{header|Ada}}==
[[Ada]] has build-in complex numbers of different accuracy (see also [[Arithmetic/Complex]]. The solution uses these types. It is thus generic, parametrized by the given complex type:
<lang Ada>
with Ada.Numerics.Generic_Complex_Arrays;
with Ada.Numerics.Generic_Complex_Elementary_Functions;

generic
with package Arrays is
new Ada.Numerics.Generic_Complex_Arrays (<>);
use Arrays;
with package Elementary_Functions is
new Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Types);
function Generic_FFT (X : Complex_Vector) return Complex_Vector;
</lang>
Generic_FFT is a generic function which takes as the parameters instances of generic packages Ada.Numerics.Generic_Complex_Arrays ([http://www.adaic.org/resources/add_content/standards/05rm/html/RM-G-3-2.html ARM G.3.2]) for complex vectors and Ada.Numerics.Generic_Complex_Elementary_Functions ([http://www.adaic.org/resources/add_content/standards/05rm/html/RM-G-1-2.html ARM G.1.2]) for elementary functions like exp. The implementation is:
<lang Ada>
with Ada.Numerics; use Ada.Numerics;

function Generic_FFT (X : Complex_Vector) return Complex_Vector is
use Complex_Types;
use Elementary_Functions;
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
A : constant Complex := Even (Even'First + K);
B : constant Complex := Odd (Odd'First + K);
begin
Even (Even'First + K) := A + B / F ** K;
Odd (Odd'First + K) := A - B / F ** K;
end;
end loop;
return Even & Odd;
end;
end if;
end FFT;
begin
return FFT (X, X'Length, 1);
end Generic_FFT;
</lang>
The following is an example of use:
<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 Test_FFT is
function FFT is
new Generic_FFT
( Ada.Numerics.Complex_Arrays,
Ada.Numerics.Complex_Elementary_Functions
);
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 Test_FFT;
</lang>
Sample output:
<pre>
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)
</pre>
=={{header|D}}==
=={{header|D}}==
<lang d>import std.stdio, std.math, std.algorithm, std.range;
<lang d>import std.stdio, std.math, std.algorithm, std.range;
Line 23: Line 111:
}</lang>
}</lang>
Output:<pre>[4+0i, 1+-2.41421i, 0+0i, 1+-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]</pre>
Output:<pre>[4+0i, 1+-2.41421i, 0+0i, 1+-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]</pre>

=={{header|Go}}==
=={{header|Go}}==
<lang go>package main
<lang go>package main

Revision as of 15:10, 20 March 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

Ada has build-in complex numbers of different accuracy (see also Arithmetic/Complex. The solution uses these types. It is thus generic, parametrized by the given complex type: <lang Ada> with Ada.Numerics.Generic_Complex_Arrays; with Ada.Numerics.Generic_Complex_Elementary_Functions;

generic 

  with package Arrays is
     new Ada.Numerics.Generic_Complex_Arrays (<>);
  use Arrays;
  with package Elementary_Functions is
     new Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Types);

function Generic_FFT (X : Complex_Vector) return Complex_Vector; </lang> Generic_FFT is a generic function which takes as the parameters instances of generic packages Ada.Numerics.Generic_Complex_Arrays (ARM G.3.2) for complex vectors and Ada.Numerics.Generic_Complex_Elementary_Functions (ARM G.1.2) for elementary functions like exp. The implementation is: <lang Ada> with Ada.Numerics; use Ada.Numerics;

function Generic_FFT (X : Complex_Vector) return Complex_Vector is 

  use Complex_Types;
  use Elementary_Functions;
  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
                 A : constant Complex := Even (Even'First + K);
                 B : constant Complex := Odd  (Odd'First  + K);
              begin
                 Even (Even'First + K) := A + B / F ** K;
                 Odd  (Odd'First  + K) := A - B / F ** K;
              end;
           end loop;
           return Even & Odd;
        end;
     end if;
  end FFT;

begin 

  return FFT (X, X'Length, 1);

end Generic_FFT; </lang> The following is an example of use: <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 Test_FFT is 

  function FFT is
     new Generic_FFT
         (  Ada.Numerics.Complex_Arrays,
            Ada.Numerics.Complex_Elementary_Functions
         );
  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 Test_FFT; </lang> Sample 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)

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));
 creal[] l, r;
 foreach (k; 0 .. N/2) {
   l ~= ev[k] + expi(-2 * PI * k / N) * od[k];
   r ~= ev[k] - expi(-2 * PI * k / N) * od[k];
 }
 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]

Go

<lang go>package main

import ( 

   "fmt"
   "math"
   "cmath"

)

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 := cmath.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)

J

Based on j:Essays/FFT, with some simplifications, sacrificing accuracy, optimizations and convenience not visible here, 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.

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