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>

{{trans|Fortran}}

{{works with|Dyalog APL}}

T←even×*(0J¯2×(○1)×(¯1+⍳k)÷N)

(odd+T),odd-T

 

'''Example:'''

 

{{out}}

 

=={{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>

}

return 0;

{{out}}

<pre>

</pre>

 

 

;;; This is adapted from the Python sample; it uses lists for simplicity.

;;; Production code would use complex arrays (for compiler optimization).

;; finally, concatenate the sum and difference of the lists

(return (mapcan #'mapcar '(+ -) `(,ffta ,ffta) `(,aux ,aux)))))))

{{out}}

;;; Demonstrates printing an FFT in both rectangular and polar form:

CL-USER> (mapc (lambda (c) (format t "~&~6F~6@Fi = ~6Fe^~6@Fipi"

#C(0.9999999999999999D0 0.4142135623730949D0) #C(0.0D0 0.0D0)

#C(0.9999999999999997D0 2.414213562373095D0))

 

=={{header|Crystal}}==

{{trans|Ruby}}

 

def fft(x : Array(Int32 | Float64)) #: Array(Complex)

return [x[0].to_c] if x.size <= 1

end

fft([1,1,1,1,0,0,0,0]).each{ |c| puts c }

{{out}}

<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

 

An implementation of the radix-2 algorithm, which works for any vector for length that is a power of 2:

function fft(a)

y1 = Any[]; y2 = Any[]

return vcat(y1,y2)

end

 

=={{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}}

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

 

 

1) the function fft

A more usefull example can be seen in http://lambdaway.free.fr/lambdaspeech/?view=zorg

 

 

=={{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 halfn = n div 2

 

result[k] = even[k] + a

result[k + halfn] = even[k] - a

 

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>

 

=={{header|Perl}}==

use warnings;

use Math::Complex;

}

 

{{out}}

<pre>4

0

1+2.41421356237309i</pre>

 

=={{header|Phix}}==

{{out}}

<pre>

 

Complex Class File:

<?php

 

}

 

 

Example:

<?php

 

echo PHP_EOL;

 

 

 

=={{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}}==

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

numeric digits length( pi() ) - length(.) /*limited by the PI function result. */

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

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

pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862

{{out|output|text=&nbsp; when using the default inputs of: &nbsp; &nbsp; <tt> 1 &nbsp; 1 &nbsp; 1 &nbsp; 1 &nbsp; 0 </tt>}}

<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|Rust}}==

{{trans|C}}

use num::complex::Complex;

use std::f64::consts::PI;

let output = fft(&input);

show("output:", &output);

{{out}}

<pre>

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

</pre>

=={{header|Scala}}==

{{libheader|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|SequenceL}}==

import <Utilities/Math.sl>;

import <Utilities/Sequence.sl>;

x when n <= 1

else

 

{{out}}

cmd:>fft(makeComplex([1,1,1,1,0,0,0,0],0))

[(Imaginary:0.00000000,Real:4.00000000),(Imaginary:-2.41421356,Real:1.00000000),(Imaginary:0.00000000,Real:0.00000000),(Imaginary:-0.41421356,Real:1.00000000),(Imaginary:0.00000000,Real:0.00000000),(Imaginary:0.41421356,Real:1.00000000),(Imaginary:0.00000000,Real:0.00000000),(Imaginary:2.41421356,Real:1.00000000)]

</pre>

 

See the '''[https://www.stata.com/help.cgi?mf_fft fft function]''' in Mata help, and in the FAQ: [https://www.stata.com/support/faqs/mata/discrete-fast-fourier-transform/ How can I calculate the Fourier coefficients of a discretely sampled function in Stata?].

 

: a=1,2,3,4

: fft(a)

1 | 10 -2 - 2i -2 -2 + 2i |

+-----------------------------------------+

 

=== fft command ===

Stata can also compute FFT using the undocumented '''fft''' command. Here is an example showing its syntax. A time variable must have been set prior to calling this command. Notice that in order to get the same result as Mata's fft() function, in both the input and the output variables the imaginary part must be passed '''first'''.

 

set obs 4

gen t=_n

tsset t

fft y x, gen(v u)

 

'''Output'''

| -2 2 |

+-----------------+</pre>

 

=={{header|SystemVerilog}}==

I could have written a more beautiful code by using non-blocking assignments in the bit_reverse_order function, but it could not be coded in a function, so FFT could not be implemented as a function as well.

 

 

package math_pkg;

 

endclass

 

Now let's perform the standard test

/// @Author: Alexandre Felipe (o.alexandre.felipe@gmail.com)

/// @Date: 2018-Jan-25

end

endmodule

By running the sanity test it outputs the following

<pre>

A more reliable test is to implement the Discrete Fourier Transform by its definition and compare the results obtained by FFT and by definition evaluation. For that let's create a class with a random data vector, and each time the vector is randomized the FFT is calculated and the output is compared by the result obtained by the definition.

 

/// @Author: Alexandre Felipe (o.alexandre.felipe@gmail.com)

/// @Date: 2018-Jan-25

endfunction

endclass

 

Now let's create a code that tests the FFT with random inputs for different sizes.

Uses a generate block since the number of samples is a parameter and must be defined at compile time.

/// @Author: Alexandre Felipe (o.alexandre.felipe@gmail.com)

/// @Date: 2018-Jan-25

endgenerate

endmodule

 

Simulating the fft_test_by_definition we get the following output:

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