Fast Fourier transform: Difference between revisions

and results in a sequence of equal length, again of complex numbers.

If you need to restrict yourself to real numbers, the output should

 

The classic version is the recursive Cooley–Tukey FFT. [http://en.wikipedia.org/wiki/Cooley–Tukey_FFT_algorithm Wikipedia] has pseudo-code for that.

Further optimizations are possible but not required.

<br><br>

 

=={{header|Ada}}==

a user instance of Ada.Numerics.Generic_Complex_Arrays.

 

with Ada.Numerics.Generic_Complex_Arrays;

use Complex_Arrays;

function Generic_FFT (X : Complex_Vector) return Complex_Vector;

 

with Ada.Numerics;

with Ada.Numerics.Generic_Complex_Elementary_Functions;

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

end Generic_FFT;

 

Example:

 

with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays;

with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;

end loop;

end;

 

{{out}}

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-2.3.5 algol68g-2.3.5].}}

{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}}

 

OP DICE = ([]SCALAR in, INT step)[]SCALAR: (

coef

FI

# -*- coding: utf-8 -*- #

 

$"1¼ cycle wave: "$, compl array fmt, one and a quarter wave ft, $l$

))

{{out}}

<pre>

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,

</pre>

 

=={{header|BBC BASIC}}==

{{works with|BBC BASIC for Windows}}

DIM Complex{r#, i#}

FOR I% = 0 TO 7

READ in{(I%)}.r#

PRINT in{(I%)}.r# "," in{(I%)}.i#

NEXT

c.i# = i#

ENDPROC

{{out}}

<pre>Input (real, imag):

Note: array size is assumed to be power of 2 and not checked by code;

you can just pad it with 0 otherwise.

 

#include <stdio.h>

}

 

{{out}}

<pre>Data: 1 1 1 1 0 0 0 0

===OS X / iOS===

OS X 10.7+ / iOS 4+

#include <Accelerate/Accelerate.h>

 

return 0;

{{out}}

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

 

=={{header|C++}}==

#include <iostream>

#include <valarray>

// Better optimized but less intuitive

// !!! Warning : in some cases this code make result different from not optimased version above (need to fix bug)

void fft(CArray &x)

{

unsigned int N = x.size(), k = N, n;

double thetaT = 3.14159265358979323846264338328L / N;

while (k > 1)

{

}

return 0;

{{out}}

<pre>

</pre>

 

 

{{out}}

 

 

 

{{out}}

<pre>

</pre>

 

=={{header|D}}==

===Standard Version===

import std.stdio, std.numeric;

 

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

{{out}}

<pre>[4+0i, 1-2.41421i, 0-0i, 1-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]</pre>

===creals Version===

Built-in complex numbers will be deprecated.

 

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

void main() @safe {

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

{{out}}

<pre>[4+0i, 1+-2.41421i, 0+0i, 1+-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]</pre>

 

===Phobos Complex Version===

 

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

void main() {

[1.0, 1, 1, 1, 0, 0, 0, 0].map!complex.array.fft.writeln;

{{out}}

<pre>[4+0i, 1-2.41421i, 0+0i, 1-0.414214i, 0+0i, 1+0.414214i, 0+0i, 1+2.41421i]</pre>

 

=={{header|EchoLisp}}==

(define -∏*2 (complex 0 (* -2 PI)))

 

→ #( 4+0i 1-2.414213562373095i 0+0i 1-0.4142135623730949i

0+0i 1+0.4142135623730949i 0+0i 1+2.414213562373095i)

 

=={{header|ERRE}}==

PROGRAM FFT

 

END FOR

END PROGRAM

{{out}}

<pre>

 

=={{header|Factor}}==

IN: USE math.transforms.fft

IN: { 1 1 1 1 0 0 0 0 } fft .

C{ 0.9999999999999997 2.414213562373095 }

}

 

=={{header|Fortran}}==

module fft_mod

implicit none

end do

 

{{out}}

<pre>

( 0.000000000000000, 0.000000000000000i )

( 1.000000000000000, 2.414213562373095i )</pre>

 

=={{header|GAP}}==

# In GAP, E(n) = exp(2*i*pi/n), a primitive root of the unity.

 

 

InverseFourier(last);

 

=={{header|Go}}==

 

import (

fmt.Printf("%8.4f\n", c)

}

{{out}}

<pre>

( 0.0000 +0.0000i)

( 1.0000 +2.4142i)

</pre>

 

=={{header|Haskell}}==

 

-- Cooley-Tukey

exp' k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)

 

{{out}}

 

=={{header|Idris}}==

 

import Data.Complex

 

main : IO()

 

{{out}}

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

 

rou =: ^@j.@o.@(% #)@i.@-: NB. roots of unity

floop =: 4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.'

 

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

 

0 0.92388 0.707107 0.382683 1 0.382683 0.707107 0.92388 0 0.92388 0.707107 0.382683 1 0.382683 0.707107 0.92388

 

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

Im=: {:@+.@fft

M=: 4#1 0

4 1 0 1 0 1 0 1

Im M

 

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

=={{header|Java}}==

{{trans|C sharp}}

 

public class FastFourierTransform {

return String.format("(%f,%f)", re, im);

}

 

<pre>Results:

variants on this page.

 

complex fast fourier transform and inverse from

http://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B

//test code

//console.log( cfft([1,1,1,1,0,0,0,0]) );

Very very basic Complex number that provides only the components

required by the code above.

basic complex number arithmetic from

http://rosettacode.org/wiki/Fast_Fourier_transform#Scala

else

console.log(this.re.toString()+'+'+this.im.toString()+'j');

 

=={{header|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):

 

# e(ix) = cos(x) + i sin(x)

====FFT====

length as $N

| if $N <= 1 then .

| [ range(0; $N/2) | cplus($even[.]; cmult( expi(-2*$pi*./$N); $odd[.] )) ] +

[ range(0; $N/2) | cminus($even[.]; cmult( expi(-2*$pi*./$N); $odd[.] )) ]

Example:

{{Out}}

[[4,-0],[1,-2.414213562373095],

 

=={{header|Julia}}==

{{out}}

4.0+0.0im

1.0-2.41421im

1.0+0.414214im

0.0+0.0im

 

=={{header|Kotlin}}==

From Scala.

 

class Complex(val re: Double, val im: Double) {

private val a = "%1.3f".format(re)

private val b = "%1.3f".format(abs(im))

 

fun fft(a: Array<Complex>) = _fft(a, Complex(0.0, 2.0), 1.0)

fun rfft(a: Array<Complex>) = _fft(a, Complex(0.0, -2.0), 2.0)

left + right

}

 

 

fun main(args: Array<String>) {

a.println()

FFT.rfft(a).println()

 

{{out}}

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

[1.000, 1.000, 1.000, 1.000, 0.000, 2.000i, 0.000, 0.000]</pre>

 

=={{header|Liberty BASIC}}==

P =8

S =int( log( P) /log( 2) +0.9999)

 

end

M Re( M) Im( M)

 

=={{header|Lua}}==

complex = {__mt={} }

data = toComplex{1, 1, 1, 1, 0, 0, 0, 0};

 

 

=={{header|Maple}}==

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

 

with( DiscreteTransforms ):

 

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

 

<pre>

</pre>

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

v := Vector( [1,1,1,1,0,0,0,0], datatype=complex[8] ):

 

 

v;

 

<pre>

[1. + 2.41421356237309 I ]

</pre>

InverseFourierTransform( v, normalization=full, inplace ):

 

v;

 

<pre>

produce output that is consistent with most other FFT routines.

 

Fourier[{1,1,1,1,0,0,0,0}, FourierParameters->{1,-1}]

 

{{out}}

<pre>{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}</pre>

 

=={{header|MATLAB}} / {{header|Octave}}==

Matlab/Octave have a builtin FFT function.

 

{{out}}

<pre>ans =

 

=={{header|Maxima}}==

fft([1, 2, 3, 4]);

 

=={{header|Nim}}==

{{trans|Python}}

 

# Works with floats and complex numbers as input

let n = x.len

return

var evens, odds = newSeq[T]()

for i, v in x:

var (even, odd) = (fft(evens), fft(odds))

 

 

 

for i in fft(@[1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]):

{{out}}

<pre>4.000

2.613

1.082

1.082

2.613</pre>

 

=={{header|OCaml}}==

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

 

let fac k n =

let indata = [one;one;one;one;zero;zero;zero;zero] in

show indata;

 

{{out}}

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

</pre>

 

=={{header|ooRexx}}==

{{trans|PL/I}} Output as shown in REXX

list='1 1 1 1 0 0 0 0'

n=words(list)

else return "-" value~abs

 

{{out}}

<pre>---data--- num real-part imaginary-part

=={{header|PARI/GP}}==

Naive implementation, using the same testcase as Ada:

{{out}}

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

{{out}}

<pre>

{ 0.000000, 0.000000},

{ 0.000000, 0.000000}}

</pre>

 

=={{header|PicoLisp}}==

{{works with|PicoLisp|3.1.0.3}}

 

(scl 4)

(native "libfftw3.so" "fftw_destroy_plan" NIL P)

(native "libfftw3.so" "fftw_free" NIL Out)

Test:

(tab (6 8)

(round (car R))

{{out}}

<pre> 4.000 0.000

 

=={{header|PL/I}}==

 

/* In-place Cooley-Tukey FFT */

end;

 

{{out}}

<pre> 4.000000000000+0.000000000000I

0.000000000000+0.000000000000I

0.999999999999+2.414213562373I</pre>

 

=={{header|Prolog}}==

{{trans|Python}}Note: Similar algorithmically to the python example.

{{works with|SWI Prolog|Version 6.2.6 by Jan Wielemaker, University of Amsterdam}}

%_______________________________________________________________

% Arithemetic for complex numbers; only the needed rules

write(R), (I>=0, write('+'),fail;write(I)), write('j, '),

fail; write(']'), nl).

 

{{out}}

=={{header|Python}}==

===Python: Recursive===

 

def fft(x):

 

print( ' '.join("%5.3f" % abs(f)

 

{{out}}

 

===Python: Using module [http://numpy.scipy.org/ numpy]===

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

 

=={{header|R}}==

The function "fft" is readily available in R

{{out}}

<pre>4+0.000000i 1-2.414214i 0+0.000000i 1-0.414214i 0+0.000000i 1+0.414214i 0+0.000000i 1+2.414214i</pre>

 

=={{header|Racket}}==

#lang racket

(require math)

(array-fft (array #[1. 1. 1. 1. 0. 0. 0. 0.]))

 

{{out}}

0.9999999999999997+2.414213562373095i])

</pre>

 

=={{header|REXX}}==

<br>('''cos'''ine &nbsp; and &nbsp; reduce radians to a unit circle).

 

 

This REXX program also adds zero values &nbsp; if &nbsp; the number of data points in the list doesn't exactly equal to a

<br>power of two. &nbsp; This is known as &nbsp; ''zero-padding''.

arg data /*ARG verb uppercases the DATA from CL.*/

/* [↓] TRANSLATE allows I & J*/ /*║ Numbers in data can be in ║*/

do j=0 for size /*║ seven formats: real ║*/

parse var _ #.1.j '' $ 1 "," #.2.j /*║ ,imag ║*/

if $=='J' then parse var #.1.j #2.j "J" #.1.j /*║ nnnJ ║*/

end /*j*/

say

say

do p-ph; halfPtr=ptr % 2

end /*p-ph*/

 

end /*j*/

end /*j*/ /* [↑] swap two sets of values. */

end /*i*/

end /*j*/

end /*k*/

end /*p*/

call hdr

exit /*stick a fork in it, we're all done. */

/*──────────────────────────────────────────────────────────────────────────────────────*/

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

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862

<pre>

 

 

</pre>

 

=={{header|Ruby}}==

return vec if vec.size <= 1

evens_odds = vec.partition.with_index{|_,i| i.even?}

end

{{Out}}

<pre>

 

=={{header|Run BASIC}}==

sig = int(log(cnt) /log(2) +0.9999)

 

print " "; i;" ";using("##.#",rel(i));" ";img(i)

next i

<pre> Num real imag

0 4.0 0

6 0.0 0

7 1.0 2.41421356</pre>

 

=={{header|Scala}}==

{{works with|Scala|2.10.4}}

Imports and Complex arithmetic:

 

case class Complex(re: Double, im: Double) {

val r = (cosh(c.re) + sinh(c.re))

Complex(cos(c.im), sin(c.im)) * r

 

The FFT definition itself:

if (cSeq.length == 1) {

return cSeq

 

def fft(cSeq: Seq[Complex]): Seq[Complex] = _fft(cSeq, Complex(0, 2), 1)

 

Usage:

Complex(0,0), Complex(0,2), Complex(0,0), Complex(0,0))

 

println(fft(data))

 

{{out}}

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

 

=={{header|Scilab}}==

Scilab has a builtin FFT function.

 

 

=={{header|Sidef}}==

{{trans|Perl}}

arr.len == 1 && return arr

 

var wave = sequence.map {|n| ::sin(n * Num.tau / sequence.len * cycles) }

say "wave:#{wave.map{|w| '%6.3f' % w }.join(' ')}"

{{out}}

<pre>

</pre>

 

 

 

{{out}}

</pre>

 

=={{header|Tcl}}==

{{tcllib|math::constants}}

{{tcllib|math::fourier}}

package require math::fourier

 

# Convert to magnitudes for printing

set fft2 [waveMagnitude $fft]

{{out}}

<pre>

=={{header|Ursala}}==

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

#import flo

 

#cast %jLW

 

{{out}}

<pre>(

0.000e+00+0.000e+00j,

1.000e+00+2.414e+00j>)</pre>

 

=={{header|zkl}}==

v:=GSL.ZVector(8).set(1,1,1,1);

{{out}}

<pre>