Fast Fourier transform: Difference between revisions

{{trans|Python}}

 

V n = x.len

I n <= 1

(0 .< n I/ 2).map(k -> @even[k] - @t[k])

 

 

{{out}}

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

1>k←2÷⍨N←⍴⍵:⍵

0≠1|2⍟N:'Argument must be a power of 2 in length'

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

 

=={{header|C sharp|C#}}==

using System.Numerics;

using System.Linq;

}

}

{{out}}

<pre>

 

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

#include <iostream>

#include <valarray>

}

return 0;

{{out}}

<pre>

and condenses also the inverse part, by a keyword.

 

(defun fft (a &key (inverse nil) &aux (n (length a)))

"Perform the FFT recursively on input vector A.

:do (setq ⍵ (* ⍵ ⍵_n))

:finally (return â)))))))

From here on the old solution.

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

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>

{{libheader| System.Math}}

{{Trans|C#}}

program Fast_Fourier_transform;

 

writeln(c);

readln;

{{out}}

<pre>4 + 0i

 

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

 

=={{header|FreeBASIC}}==

'press any key for the next image.

'131072 samples: the FFT is fast indeed.

cls

next i

<b>(Images only)</b>

 

=={{header|Frink}}==

Frink has a built-in FFT function that can produce results based on different conventions. The following is not the default convention, but matches many of the other results in this page.

println[joinln[format[a, 1, 5]]]

{{out}}

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

 

=={{header|Golfscript}}==

 

{.,.({[\.2%fft\(;2%fft@-1?-1\?-2?:w;.,,{w\?}%[\]zip{{*}*}%]zip.{{+}*}%\{{-}*}%+}{;}if}:fft;

 

[1 1 1 1 0 0 0 0]fft n*

{{out}}

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

f::{p::2:#x;e::ff2(*'p);o::ff2({x@1}'p);k::-1;

t::{k::k+1;cmul(cexp(cdiv(cmul([0 -2];(k*pi),0);n,0));x)}'o;

(e cadd't),e csub't};

:[n<2;x;f(x)]};

Example (rounding to 4 decimal digits):

{{out}}

<pre>[[4.0 0.0]

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

-- Beginning with Lua 5.2 you have to write

print("orig:", table.unpack(data))

 

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

Here is a user-space definition for good measure.

 

fft[l__] :=

Join[#, #] &@fft@l[[1 ;; ;; 2]] +

fft[l[[2 ;; ;; 2]]])

 

{{out}}

<pre>4.

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

 

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

{{out}}

<pre>4.000

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

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

=== Recursive ===

 

 

 

=={{header|Perl}}==

{{trans|Raku}}

use warnings;

use Math::Complex;

}

 

{{out}}

<pre>4

 

=={{header|Phix}}==

<span style="color: #000080;font-style:italic;">--

-- demo\rosetta\FastFourierTransform.exw

<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #008000;">"\nResults of %d-point inverse fft (rounded to 6 d.p.):\n\n"</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">))</span>

<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ifft</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fft</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)))</span>

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

$Len = $Arr.Count

 

Return $Output

}

{{out}}

<pre>

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

=={{header|Raku}}==

(formerly Perl 6)

return @_ if @_ == 1;

my @evn = fft( @_[0, 2 ... *] );

return flat @evn »+« @odd, @evn »-« @odd;

}

{{out}}

<pre>4+0i

0+0i

0+0i

0+0i

 

For the fun of it, here is a purely functional version:

 

@_ == 1 ?? @_ !!

fft(@_[0,2...*]) «+«

fft(@_[1,3...*]) «*« map &cis, (0,-τ/@_...^-τ)

 

This particular version is numerically inaccurate though, because of the pi approximation. It is possible to fix it with the 'cisPi' function available in the TrigPi module:

 

use TrigPi;

@_ == 1 ?? @_ !!

fft(@_[0,2...*]) «+«

fft(@_[1,3...*]) «*« map &cisPi, (0,-2/@_...^-2)

 

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

Programming note: &nbsp; the numeric precision (decimal digits) is only restricted by the number of decimal digits in the &nbsp;

<br>'''pi''' &nbsp; variable &nbsp; (which is defined in the penultimate assignment statement in the REXX program.

 

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

=={{header|Rust}}==

{{trans|C}}

use num::complex::Complex;

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

let output = fft(&input);

show("output:", &output);

{{out}}

<pre>

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

import <Utilities/Math.sl>;

import <Utilities/Sequence.sl>;

x when n <= 1

else

 

{{out}}

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

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

{{trans|Kotlin}}

 

import Numerics

 

 

print(fft(dat).map({ $0.pretty }))

 

{{out}}

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

{{libheader|Wren-complex}}

{{libheader|Wren-fmt}}

 

var ditfft2 // recursive

for (i in 0...y.count) y[i] = Complex.zero

ditfft2.call(x, y, x.count, 1)

 

{{out}}

 

=={{header|zkl}}==

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

{{out}}

<pre>