Fast Fourier transform: Difference between revisions
m
→Recursive
m (→Recursive) |
m (→Recursive) |
||
(9 intermediate revisions by 6 users not shown) | |||
Line 17:
{{trans|Python}}
<
V n = x.len
I n <= 1
Line 27:
(0 .< n I/ 2).map(k -> @even[k] - @t[k])
print(fft([Complex(1.0), 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]).map(f -> ‘#1.3’.format(abs(f))).join(‘ ’))</
{{out}}
Line 39:
a user instance of Ada.Numerics.Generic_Complex_Arrays.
<syntaxhighlight lang="ada">
with Ada.Numerics.Generic_Complex_Arrays;
Line 47:
use Complex_Arrays;
function Generic_FFT (X : Complex_Vector) return Complex_Vector;
</
<syntaxhighlight lang="ada">
with Ada.Numerics;
with Ada.Numerics.Generic_Complex_Elementary_Functions;
Line 89:
return FFT (X, X'Length, 1);
end Generic_FFT;
</syntaxhighlight>
Example:
<syntaxhighlight lang="ada">
with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays;
with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;
Line 114:
end loop;
end;
</syntaxhighlight>
{{out}}
Line 134:
{{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''.}}
'''File: Template.Fast_Fourier_transform.a68'''<
OP DICE = ([]SCALAR in, INT step)[]SCALAR: (
Line 171:
coef
FI
);</
# -*- coding: utf-8 -*- #
Line 192:
$"1¼ cycle wave: "$, compl array fmt, one and a quarter wave ft, $l$
))
)</
{{out}}
<pre>
Line 202:
{{trans|Fortran}}
{{works with|Dyalog APL}}
<
1>k←2÷⍨N←⍴⍵:⍵
0≠1|2⍟N:'Argument must be a power of 2 in length'
Line 209:
T←even×*(0J¯2×(○1)×(¯1+⍳k)÷N)
(odd+T),odd-T
}</
'''Example:'''
<
{{out}}
Line 220:
=={{header|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
<
DIM Complex{r#, i#}
Line 229:
FOR I% = 0 TO 7
READ in{(I%)}.r#
out{(I%)}.r# = in{(I%)}.r#
PRINT in{(I%)}.r# "," in{(I%)}.i#
NEXT
Line 264 ⟶ 265:
c.i# = i#
ENDPROC
</syntaxhighlight>
{{out}}
<pre>Input (real, imag):
Line 289 ⟶ 290:
Note: array size is assumed to be power of 2 and not checked by code;
you can just pad it with 0 otherwise.
Also, <code>complex</code> is C99 standard.<
#include <stdio.h>
Line 342 ⟶ 343:
}
</syntaxhighlight>
{{out}}
<pre>Data: 1 1 1 1 0 0 0 0
Line 349 ⟶ 350:
===OS X / iOS===
OS X 10.7+ / iOS 4+
<
#include <Accelerate/Accelerate.h>
Line 389 ⟶ 390:
return 0;
}</
{{out}}
<pre>Data: 1 1 1 1 0 0 0 0
Line 395 ⟶ 396:
=={{header|C sharp|C#}}==
<
using System.Numerics;
using System.Linq;
Line 497 ⟶ 498:
}
}
}</
{{out}}
<pre>
Line 512 ⟶ 513:
=={{header|C++}}==
<
#include <iostream>
#include <valarray>
Line 636 ⟶ 637:
}
return 0;
}</
{{out}}
<pre>
Line 668 ⟶ 669:
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.
Line 707 ⟶ 708:
:do (setq ⍵ (* ⍵ ⍵_n))
:finally (return â)))))))
</syntaxhighlight>
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).
Line 725 ⟶ 726:
;; finally, concatenate the sum and difference of the lists
(return (mapcan #'mapcar '(+ -) `(,ffta ,ffta) `(,aux ,aux)))))))
</syntaxhighlight>
{{out}}
<
;;; Demonstrates printing an FFT in both rectangular and polar form:
CL-USER> (mapc (lambda (c) (format t "~&~6F~6@Fi = ~6Fe^~6@Fipi"
Line 746 ⟶ 747:
#C(0.9999999999999999D0 0.4142135623730949D0) #C(0.0D0 0.0D0)
#C(0.9999999999999997D0 2.414213562373095D0))
</syntaxhighlight>
=={{header|Crystal}}==
{{trans|Ruby}}
<
def fft(x : Array(Int32 | Float64)) #: Array(Complex)
Line 762 ⟶ 763:
fft([1,1,1,1,0,0,0,0]).each{ |c| puts c }
</syntaxhighlight>
{{out}}
<pre>
Line 777 ⟶ 778:
=={{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>
Line 787 ⟶ 788:
===creals Version===
Built-in complex numbers will be deprecated.
<
const(creal)[] fft(in creal[] x) pure /*nothrow*/ @safe {
Line 801 ⟶ 802:
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*/ {
Line 821 ⟶ 822:
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>
Line 829 ⟶ 830:
{{libheader| System.Math}}
{{Trans|C#}}
<syntaxhighlight lang="delphi">
program Fast_Fourier_transform;
Line 917 ⟶ 918:
writeln(c);
readln;
end.</
{{out}}
<pre>4 + 0i
Line 929 ⟶ 930:
=={{header|EchoLisp}}==
<
(define -∏*2 (complex 0 (* -2 PI)))
Line 947 ⟶ 948:
→ #( 4+0i 1-2.414213562373095i 0+0i 1-0.4142135623730949i
0+0i 1+0.4142135623730949i 0+0i 1+2.414213562373095i)
</syntaxhighlight>
=={{header|ERRE}}==
<syntaxhighlight lang="erre">
PROGRAM FFT
Line 1,050 ⟶ 1,051:
END FOR
END PROGRAM
</syntaxhighlight>
{{out}}
<pre>
Line 1,069 ⟶ 1,070:
=={{header|Factor}}==
<syntaxhighlight lang="factor">
IN: USE math.transforms.fft
IN: { 1 1 1 1 0 0 0 0 } fft .
Line 1,082 ⟶ 1,083:
C{ 0.9999999999999997 2.414213562373095 }
}
</syntaxhighlight>
=={{header|Fortran}}==
<syntaxhighlight lang="fortran">
module fft_mod
implicit none
Line 1,143 ⟶ 1,144:
end do
end program test</
{{out}}
<pre>
Line 1,156 ⟶ 1,157:
=={{header|FreeBASIC}}==
<
'press any key for the next image.
'131072 samples: the FFT is fast indeed.
Line 1,374 ⟶ 1,375:
cls
next i
end</
<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]]]
</syntaxhighlight>
{{out}}
<pre>
Line 1,395 ⟶ 1,396:
=={{header|GAP}}==
<
# In GAP, E(n) = exp(2*i*pi/n), a primitive root of the unity.
Line 1,417 ⟶ 1,418:
InverseFourier(last);
# [ 1, 1, 1, 1, 0, 0, 0, 0 ]</
=={{header|Go}}==
<
import (
Line 1,448 ⟶ 1,449:
fmt.Printf("%8.4f\n", c)
}
}</
{{out}}
<pre>
Line 1,462 ⟶ 1,463:
=={{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*
</syntaxhighlight>
{{out}}
<pre>
Line 1,481 ⟶ 1,482:
=={{header|Haskell}}==
<
-- Cooley-Tukey
Line 1,497 ⟶ 1,498:
exp' k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)
main = mapM_ print $ fft [1,1,1,1,0,0,0,0]</
{{out}}
Line 1,511 ⟶ 1,512:
=={{header|Idris}}==
<syntaxhighlight lang="text">module Main
import Data.Complex
Line 1,541 ⟶ 1,542:
main : IO()
main = traverse_ printLn $ fft [1,1,1,1,0,0,0,0]</
{{out}}
Line 1,558 ⟶ 1,559:
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.'
fft =: ] floop&.cube rou@#</
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
0 0 0 0j8 0 0 0 0 0 0 0 0 0 0j8 0 0</
Here is a representation of an example which appears in some of the other implementations, here:
<
Im=: {:@+.@fft
M=: 4#1 0
Line 1,578 ⟶ 1,579:
4 1 0 1 0 1 0 1
Im M
0 2.41421 0 0.414214 0 _0.414214 0 _2.41421</
Note that Re and Im are not functions of 1 and 0 but are functions of the complete sequence.
Line 1,586 ⟶ 1,587:
=={{header|Java}}==
{{trans|C sharp}}
<
public class FastFourierTransform {
Line 1,680 ⟶ 1,681:
return String.format("(%f,%f)", re, im);
}
}</
<pre>Results:
Line 1,696 ⟶ 1,697:
variants on this page.
<
complex fast fourier transform and inverse from
http://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B
Line 1,761 ⟶ 1,762:
//test code
//console.log( cfft([1,1,1,1,0,0,0,0]) );
//console.log( icfft(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
Line 1,813 ⟶ 1,814:
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====
<syntaxhighlight lang="jq">
# multiplication of real or complex numbers
def cmult(x; y):
Line 1,841 ⟶ 1,842:
# e(ix) = cos(x) + i sin(x)
def expi(x): [ (x|cos), (x|sin) ];</
====FFT====
<
length as $N
| if $N <= 1 then .
Line 1,851 ⟶ 1,852:
| [ range(0; $N/2) | cplus($even[.]; cmult( expi(-2*$pi*./$N); $odd[.] )) ] +
[ range(0; $N/2) | cminus($even[.]; cmult( expi(-2*$pi*./$N); $odd[.] )) ]
end;</
Example:
<
</syntaxhighlight>
{{Out}}
[[4,-0],[1,-2.414213562373095],
Line 1,862 ⟶ 1,863:
=={{header|Julia}}==
<
fft([1,1,1,1,0,0,0,0])</
{{out}}
<
4.0+0.0im
1.0-2.41421im
Line 1,874 ⟶ 1,875:
1.0+0.414214im
0.0+0.0im
1.0+2.41421im</
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[]
Line 1,893 ⟶ 1,894:
return vcat(y1,y2)
end
</syntaxhighlight>
=={{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)]};
n::#x;k::{(2^x)<n}{1+x}:~1;n#ff2({x,0}'x,&(2^k)-n)}</
Example (rounding to 4 decimal digits):
<
{{out}}
<pre>[[4.0 0.0]
Line 1,916 ⟶ 1,917:
=={{header|Kotlin}}==
From Scala.
<
class Complex(val re: Double, val im: Double) {
Line 1,935 ⟶ 1,936:
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)
Line 1,964 ⟶ 1,965:
left + right
}
}</
<
fun main(args: Array<String>) {
Line 1,975 ⟶ 1,976:
a.println()
FFT.rfft(a).println()
}</
{{out}}
Line 1,982 ⟶ 1,983:
=={{header|Lambdatalk}}==
<
1) the function fft
Line 2,093 ⟶ 2,094:
A more usefull example can be seen in http://lambdaway.free.fr/lambdaspeech/?view=zorg
</syntaxhighlight>
=={{header|Liberty BASIC}}==
<syntaxhighlight lang="lb">
P =8
S =int( log( P) /log( 2) +0.9999)
Line 2,202 ⟶ 2,203:
end
</syntaxhighlight>
M Re( M) Im( M)
Line 2,215 ⟶ 2,216:
=={{header|Lua}}==
<
complex = {__mt={} }
Line 2,286 ⟶ 2,287:
-- Beginning with Lua 5.2 you have to write
print("orig:", table.unpack(data))
print("fft:", table.unpack(fft(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.
<syntaxhighlight lang="maple">
with( DiscreteTransforms ):
FourierTransform( <1,1,1,1,0,0,0,0>, normalization=none );
</syntaxhighlight>
<pre>
Line 2,315 ⟶ 2,316:
</pre>
Optionally, the FFT may be performed inplace on a Vector of hardware double-precision complex floats.
<syntaxhighlight lang="maple">
v := Vector( [1,1,1,1,0,0,0,0], datatype=complex[8] ):
Line 2,321 ⟶ 2,322:
v;
</syntaxhighlight>
<pre>
Line 2,340 ⟶ 2,341:
[1. + 2.41421356237309 I ]
</pre>
<syntaxhighlight lang="maple">
InverseFourierTransform( v, normalization=full, inplace ):
v;
</syntaxhighlight>
<pre>
Line 2,371 ⟶ 2,372:
produce output that is consistent with most other FFT routines.
<syntaxhighlight lang="mathematica">
Fourier[{1,1,1,1,0,0,0,0}, FourierParameters->{1,-1}]
</syntaxhighlight>
{{out}}
Line 2,380 ⟶ 2,381:
Here is a user-space definition for good measure.
<
fft[l__] :=
Join[#, #] &@fft@l[[1 ;; ;; 2]] +
Line 2,386 ⟶ 2,387:
fft[l[[2 ;; ;; 2]]])
fft[{1, 1, 1, 1, 0, 0, 0, 0}] // Column</
{{out}}
<pre>4.
Line 2,402 ⟶ 2,403:
Matlab/Octave have a builtin FFT function.
<
</syntaxhighlight>
{{out}}
<pre>ans =
Line 2,417 ⟶ 2,418:
=={{header|Maxima}}==
<
fft([1, 2, 3, 4]);
[2.5, -0.5 * %i - 0.5, -0.5, 0.5 * %i - 0.5]</
=={{header|Nim}}==
{{trans|Python}}
<
# Works with floats and complex numbers as input
Line 2,450 ⟶ 2,451:
for i in fft(@[1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]):
echo formatFloat(abs(i), ffDecimal, 3)</
{{out}}
<pre>4.000
Line 2,463 ⟶ 2,464:
=={{header|OCaml}}==
This is a simple implementation of the Cooley-Tukey pseudo-code
<
let fac k n =
Line 2,487 ⟶ 2,488:
let indata = [one;one;one;one;zero;zero;zero;zero] in
show indata;
show (fft indata)</
{{out}}
Line 2,497 ⟶ 2,498:
=={{header|ooRexx}}==
{{trans|PL/I}} Output as shown in REXX
<
list='1 1 1 1 0 0 0 0'
n=words(list)
Line 2,614 ⟶ 2,615:
else return "-" value~abs
::requires rxMath library</
{{out}}
<pre>---data--- num real-part imaginary-part
Line 2,639 ⟶ 2,640:
=={{header|PARI/GP}}==
Naive implementation, using the same testcase as Ada:
<
FFT([1,1,1,1,0,0,0,0])</
{{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>
differently, and even with "graphics"
<
install( FFT, GG );
k = 7; N = 2 ^ k;
Line 2,654 ⟶ 2,655:
\\plot( i = 1, N, v[ floor(i) ] );
print("Spectrum");
plot( i = 1, N / 2 , abs( w[floor(i)] ) * 2 / N );</
{{out}}
<pre>Spectrum
Line 2,684 ⟶ 2,685:
=== Recursive ===
{{works with|Free Pascal|3.2.0 }}
<
PROGRAM RDFT;
Line 2,708 ⟶ 2,709:
ucomplex;
{$WARN 6058 off : Call to subroutine "$1" marked as inline is not inlined}
(*) Use for variants and ucomplex (*)
Line 2,872 ⟶ 2,875:
</syntaxhighlight>
JPD 2021/12/26
=={{header|Perl}}==
{{trans|Raku}}
<
use warnings;
use Math::Complex;
Line 2,892 ⟶ 2,895:
}
print "$_\n" for fft qw(1 1 1 1 0 0 0 0);</
{{out}}
<pre>4
Line 2,904 ⟶ 2,907:
=={{header|Phix}}==
<!--<
<span style="color: #000080;font-style:italic;">--
-- demo\rosetta\FastFourierTransform.exw
Line 3,035 ⟶ 3,038:
<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>
Line 3,065 ⟶ 3,068:
Complex Class File:
<syntaxhighlight lang="php">
<?php
Line 3,107 ⟶ 3,110:
}
</
Example:
<syntaxhighlight lang="php">
<?php
Line 3,213 ⟶ 3,216:
echo PHP_EOL;
</syntaxhighlight>
Line 3,254 ⟶ 3,257:
=={{header|PicoLisp}}==
{{works with|PicoLisp|3.1.0.3}}
<
(scl 4)
Line 3,273 ⟶ 3,276:
(native "libfftw3.so" "fftw_destroy_plan" NIL P)
(native "libfftw3.so" "fftw_free" NIL Out)
(native "libfftw3.so" "fftw_free" NIL In) ) ) )</
Test:
<
(tab (6 8)
(round (car R))
(round (cadr R)) ) )</
{{out}}
<pre> 4.000 0.000
Line 3,290 ⟶ 3,293:
=={{header|PL/I}}==
<
/* In-place Cooley-Tukey FFT */
Line 3,338 ⟶ 3,341:
end;
END test;</
{{out}}
<pre> 4.000000000000+0.000000000000I
Line 3,350 ⟶ 3,353:
=={{header|POV-Ray}}==
<syntaxhighlight lang="pov-ray">
//cmd: +w0 +h0 -F -D
//Stockham algorithm
Line 3,448 ⟶ 3,451:
CdebugArr(cdata)
#debug "\n"
</syntaxhighlight>
{{out}}
Line 3,495 ⟶ 3,498:
=={{header|PowerShell}}==
<
$Len = $Arr.Count
Line 3,527 ⟶ 3,530:
Return $Output
}
</syntaxhighlight>
{{out}}
<pre>
Line 3,546 ⟶ 3,549:
{{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
Line 3,583 ⟶ 3,586:
write(R), (I>=0, write('+'),fail;write(I)), write('j, '),
fail; write(']'), nl).
</syntaxhighlight>
{{out}}
Line 3,593 ⟶ 3,596:
=={{header|Python}}==
===Python: Recursive===
<
def fft(x):
Line 3,605 ⟶ 3,608:
print( ' '.join("%5.3f" % abs(f)
for f in fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])) )</
{{out}}
Line 3,611 ⟶ 3,614:
===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)) )
4.000 2.613 0.000 1.082 0.000 1.082 0.000 2.613</
=={{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.]))
</syntaxhighlight>
{{out}}
Line 3,645 ⟶ 3,648:
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|rakudo
{{improve|raku|Not numerically accurate}}
<syntaxhighlight lang="raku" line>sub fft {
return @_ if @_ == 1;
my @evn = fft( @_[0, 2 ... *] );
Line 3,653 ⟶ 3,657:
return flat @evn »+« @odd, @evn »-« @odd;
}
.say for fft <1 1 1 1 0 0 0 0>;</
{{out}}
<pre>4+0i
1-2.414213562373095i
0+0i
1-0.4142135623730949i
0+0i
0.9999999999999999+0.4142135623730949i
0+0i
0.9999999999999997+2.414213562373095i
</pre>
For the fun of it, here is a purely functional version:
<syntaxhighlight lang="raku"
@_ == 1 ?? @_ !!
fft(@_[0,2...*]) «+«
fft(@_[1,3...*]) «*« map &cis, (0,-τ/@_...^-τ)
}</
<!-- Not really helping now
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:
<syntaxhighlight lang="raku"
use TrigPi;
@_ == 1 ?? @_ !!
fft(@_[0,2...*]) «+«
fft(@_[1,3...*]) «*« map &cisPi, (0,-2/@_...^-2)
}</
-->
=={{header|REXX}}==
Line 3,694 ⟶ 3,701:
This REXX program also adds zero values if the number of data points in the list doesn't exactly equal to a
<br>power of two. This is known as ''zero-padding''.
<
numeric digits length( pi() ) - length(.) /*limited by the PI function result. */
arg data /*ARG verb uppercases the DATA from CL.*/
Line 3,761 ⟶ 3,768:
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
r2r: return arg(1) // ( pi() * 2 ) /*reduce the radians to a unit circle. */</
Programming note: the numeric precision (decimal digits) is only restricted by the number of decimal digits in the
<br>'''pi''' variable (which is defined in the penultimate assignment statement in the REXX program.
Line 3,793 ⟶ 3,800:
=={{header|Ruby}}==
<
return vec if vec.size <= 1
evens_odds = vec.partition.with_index{|_,i| i.even?}
Line 3,802 ⟶ 3,809:
end
fft([1,1,1,1,0,0,0,0]).each{|c| puts "%9.6f %+9.6fi" % c.rect}</
{{Out}}
<pre>
Line 3,816 ⟶ 3,823:
=={{header|Run BASIC}}==
<
sig = int(log(cnt) /log(2) +0.9999)
Line 3,915 ⟶ 3,922:
print " "; i;" ";using("##.#",rel(i));" ";img(i)
next i
end</
<pre> Num real imag
0 4.0 0
Line 3,928 ⟶ 3,935:
=={{header|Rust}}==
{{trans|C}}
<
use num::complex::Complex;
use std::f64::consts::PI;
Line 3,988 ⟶ 3,995:
let output = fft(&input);
show("output:", &output);
}</
{{out}}
<pre>
Line 4,002 ⟶ 4,009:
{{works with|Scala|2.10.4}}
Imports and Complex arithmetic:
<
case class Complex(re: Double, im: Double) {
Line 4,026 ⟶ 4,033:
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
Line 4,053 ⟶ 4,060:
def fft(cSeq: Seq[Complex]): Seq[Complex] = _fft(cSeq, Complex(0, 2), 1)
def rfft(cSeq: Seq[Complex]): Seq[Complex] = _fft(cSeq, Complex(0, -2), 2)</
Usage:
<
Complex(0,0), Complex(0,2), Complex(0,0), Complex(0,0))
println(fft(data))
println(rfft(fft(data)))</
{{out}}
Line 4,067 ⟶ 4,074:
=={{header|Scheme}}==
{{works with|Chez Scheme}}
<
; Decimation-in-Time (DIT) algorithm. The input is assumed to be a vector of complex
; numbers that is a power of two in length greater than zero.
Line 4,099 ⟶ 4,106:
(dft (fft-r2dit inp)))
(printf "In: ~a~%" inp)
(printf "DFT: ~a~%" dft))</
{{out}}
<pre>
Line 4,110 ⟶ 4,117:
Scilab has a builtin FFT function.
<
=={{header|SequenceL}}==
<
import <Utilities/Math.sl>;
import <Utilities/Sequence.sl>;
Line 4,130 ⟶ 4,137:
x when n <= 1
else
complexAdd(top,z) ++ complexSubtract(top,z);</
{{out}}
Line 4,140 ⟶ 4,147:
=={{header|Sidef}}==
{{trans|Perl}}
<
arr.len == 1 && return arr
Line 4,155 ⟶ 4,162:
var wave = sequence.map {|n| ::sin(n * Num.tau / sequence.len * cycles) }
say "wave:#{wave.map{|w| '%6.3f' % w }.join(' ')}"
say "fft: #{fft(wave).map { '%6.3f' % .abs }.join(' ')}"</
{{out}}
<pre>
Line 4,167 ⟶ 4,174:
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?].
<syntaxhighlight lang="text">. mata
: a=1,2,3,4
: fft(a)
Line 4,174 ⟶ 4,181:
1 | 10 -2 - 2i -2 -2 + 2i |
+-----------------------------------------+
: end</
=== 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
Line 4,186 ⟶ 4,193:
tsset t
fft y x, gen(v u)
list u v, noobs</
'''Output'''
Line 4,205 ⟶ 4,212:
{{trans|Kotlin}}
<
import Numerics
Line 4,273 ⟶ 4,280:
print(fft(dat).map({ $0.pretty }))
print(rfft(f).map({ $0.pretty }))</
{{out}}
Line 4,286 ⟶ 4,293:
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.
<syntaxhighlight lang="systemverilog">
package math_pkg;
Line 4,374 ⟶ 4,381:
endclass
</syntaxhighlight>
Now let's perform the standard test
<syntaxhighlight lang="systemverilog">
/// @Author: Alexandre Felipe (o.alexandre.felipe@gmail.com)
/// @Date: 2018-Jan-25
Line 4,401 ⟶ 4,408:
end
endmodule
</syntaxhighlight>
By running the sanity test it outputs the following
<pre>
Line 4,427 ⟶ 4,434:
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.
<syntaxhighlight lang="systemverilog">
/// @Author: Alexandre Felipe (o.alexandre.felipe@gmail.com)
/// @Date: 2018-Jan-25
Line 4,480 ⟶ 4,487:
endfunction
endclass
</syntaxhighlight>
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.
<syntaxhighlight lang="systemverilog">
/// @Author: Alexandre Felipe (o.alexandre.felipe@gmail.com)
/// @Date: 2018-Jan-25
Line 4,499 ⟶ 4,506:
endgenerate
endmodule
</syntaxhighlight>
Simulating the fft_test_by_definition we get the following output:
Line 4,529 ⟶ 4,536:
{{tcllib|math::constants}}
{{tcllib|math::fourier}}
<
package require math::fourier
Line 4,566 ⟶ 4,573:
# Convert to magnitudes for printing
set fft2 [waveMagnitude $fft]
printwave fft2</
{{out}}
<pre>
Line 4,575 ⟶ 4,582:
=={{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
Line 4,584 ⟶ 4,591:
#cast %jLW
t = (f,g)</
{{out}}
<pre>(
Line 4,606 ⟶ 4,613:
1.000e+00+2.414e+00j>)</pre>
=={{header|
{{works with|VBA}}
{{trans|BBC_BASIC}}
Written and tested in Microsoft Visual Basic for Applications 7.1 under Office 365 Excel; but is probably useable under any recent version of VBA.
<syntaxhighlight lang="vba">Option Base 0
Private Type Complex
re As Double
im As Double
End Type
Private Function cmul(c1 As Complex, c2 As Complex) As Complex
Dim ret As Complex
ret.re = c1.re * c2.re - c1.im * c2.im
ret.im = c1.re * c2.im + c1.im * c2.re
cmul = ret
End Function
Public Sub FFT(buf() As Complex, out() As Complex, begin As Integer, step As Integer, N As Integer)
Dim i As Integer, t As Complex, c As Complex, v As Complex
If step < N Then
FFT out, buf, begin, 2 * step, N
FFT out, buf, begin + step, 2 * step, N
i = 0
While i < N
t.re = Cos(-WorksheetFunction.Pi() * i / N)
t.im = Sin(-WorksheetFunction.Pi() * i / N)
c = cmul(t, out(begin + i + step))
buf(begin + (i \ 2)).re = out(begin + i).re + c.re
buf(begin + (i \ 2)).im = out(begin + i).im + c.im
buf(begin + ((i + N) \ 2)).re = out(begin + i).re - c.re
buf(begin + ((i + N) \ 2)).im = out(begin + i).im - c.im
i = i + 2 * step
Wend
End If
End Sub
' --- test routines:
Private Sub show(r As Long, txt As String, buf() As Complex)
Dim i As Integer
r = r + 1
Cells(r, 1) = txt
For i = LBound(buf) To UBound(buf)
r = r + 1
Cells(r, 1) = buf(i).re: Cells(r, 2) = buf(i).im
Next
End Sub
Sub testFFT()
Dim buf(7) As Complex, out(7) As Complex
Dim r As Long, i As Integer
buf(0).re = 1: buf(1).re = 1: buf(2).re = 1: buf(3).re = 1
r = 0
show r, "Input (real, imag):", buf
FFT out, buf, 0, 1, 8
r = r + 1
show r, "Output (real, imag):", out
End Sub
</syntaxhighlight>
{{out}}
<pre>Input (real, imag):
1 0
1 0
1 0
1 0
0 0
0 0
0 0
0 0
Output (real, imag):
4 0
1 -2.414213562
0 0
1 -0.414213562
0 0
1 0.414213562
0 0
1 2.414213562</pre>
=={{header|V (Vlang)}}==
{{trans|Go}}
<syntaxhighlight lang="v (vlang)">import math.complex
import math
fn ditfft2(x []f64, mut y []Complex, n int, s int) {
Line 4,630 ⟶ 4,720:
println("${c:8.4f}")
}
}</
{{out}}
Line 4,653 ⟶ 4,743:
{{libheader|Wren-complex}}
{{libheader|Wren-fmt}}
<
import "./fmt" for Fmt
var ditfft2 // recursive
Line 4,679 ⟶ 4,769:
for (i in 0...y.count) y[i] = Complex.zero
ditfft2.call(x, y, x.count, 1)
for (c in y) Fmt.print("$6.4z", c)</
{{out}}
Line 4,694 ⟶ 4,784:
=={{header|zkl}}==
<
v:=GSL.ZVector(8).set(1,1,1,1);
GSL.FFT(v).toList().concat("\n").println(); // in place</
{{out}}
<pre>
|