Fast Fourier transform: Difference between revisions

From Rosetta Code
Content added Content deleted
(Ada solution added)
(→‎{{header|C++}}: Added C++ implementation)
Line 89: Line 89:
( 0.000, 0.000) ( 1.000, 2.414)
( 0.000, 0.000) ( 1.000, 2.414)
</pre>
</pre>

=={{header|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)
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 = exp(Complex(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";
}
}</lang>
Output:
<pre>
(4,0)
(1,-2.41421)
(0,0)
(1,-0.414214)
(0,0)
(1,0.414214)
(0,0)
(1,2.41421)
</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;

Revision as of 22:40, 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)

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 = exp(Complex(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";
   }

}</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));
 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=171431"