Fast Fourier transform: Difference between revisions
m (→{{header|J}}) |
|||
Line 974: | Line 974: | ||
6 0.0 0 |
6 0.0 0 |
||
7 1.0 2.41421356</pre> |
7 1.0 2.41421356</pre> |
||
=={{header|Scala}}== |
|||
</lang> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
Revision as of 11:02, 7 May 2012
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
The FFT function is defined as a generic function, instantiated upon a user instance of Ada.Numerics.Generic_Complex_Arrays.
<lang Ada> with Ada.Numerics.Generic_Complex_Arrays;
generic
with package Complex_Arrays is new Ada.Numerics.Generic_Complex_Arrays (<>); use Complex_Arrays;
function Generic_FFT (X : Complex_Vector) return Complex_Vector;
</lang>
<lang Ada> with Ada.Numerics; with Ada.Numerics.Generic_Complex_Elementary_Functions;
function Generic_FFT (X : Complex_Vector) return Complex_Vector is
package Complex_Elementary_Functions is new Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Arrays.Complex_Types);
use Ada.Numerics; use Complex_Elementary_Functions; use Complex_Arrays.Complex_Types; 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 T : constant Complex := Odd (Odd'First + K) / F ** K; begin Odd (Odd'First + K) := Even (Even'First + K) - T; Even (Even'First + K) := Even (Even'First + K) + T; end; end loop; return Even & Odd; end; end if; end FFT;
begin
return FFT (X, X'Length, 1);
end Generic_FFT; </lang>
Example:
<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 Example is
function FFT is new Generic_FFT (Ada.Numerics.Complex_Arrays); 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; </lang>
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)
ALGOL 68
Note: This specimen retains the original Python coding style.
File: Template.Fast_Fourier_transform.a68<lang algol68>PRIO DICE = 9; # ideally = 11 #
OP DICE = ([]SCALAR in, INT step)[]SCALAR: (
### Dice the array, extract array values a "step" apart ### IF step = 1 THEN in ELSE INT upb out := 0; [(UPB in-LWB in)%step+1]SCALAR out; FOR index FROM LWB in BY step TO UPB in DO out[upb out+:=1] := in[index] OD; out[@LWB in] FI
);
PROC fft = ([]SCALAR in t)[]SCALAR: (
### The Cooley-Tukey FFT algorithm ### IF LWB in t >= UPB in t THEN in t[@0] ELSE []SCALAR t = in t[@0]; INT n = UPB t + 1, half n = n % 2; [LWB t:UPB t]SCALAR coef;
[]SCALAR even = fft(t DICE 2), odd = fft(t[1:]DICE 2);
COMPL i = 0 I 1;
REAL w = 2*pi / n; FOR k FROM LWB t TO half n-1 DO COMPL cis t = scalar exp(0 I (-w * k))*odd[k]; coef[k] := even[k] + cis t; coef[k + half n] := even[k] - cis t OD; coef FI
);</lang>File: test.Fast_Fourier_transform.a68<lang algol68>#!/usr/local/bin/a68g --script #
- -*- coding: utf-8 -*- #
MODE SCALAR = COMPL; PROC (COMPL)COMPL scalar exp = complex exp; PR READ "Template.Fast_Fourier_transform.a68" PR
FORMAT real fmt := $g(0,3)$; FORMAT real array fmt := $f(real fmt)", "$; FORMAT compl fmt := $f(real fmt)"⊥"f(real fmt)$; FORMAT compl array fmt := $f(compl fmt)", "$;
test:(
[]COMPL tooth wave ft = fft((1, 1, 1, 1, 0, 0, 0, 0)), one and a quarter wave ft = fft((0, 0.924, 0.707,-0.383,-1,-0.383, 0.707, 0.924, 0,-0.924,-0.707, 0.383, 1, 0.383,-0.707,-0.924)); printf(( $"Tooth wave: "$,compl array fmt, tooth wave ft, $l$, $"1¼ cycle wave: "$, compl array fmt, one and a quarter wave ft, $l$ ))
)</lang>Output:
Tooth wave: 4.000⊥.000, 1.000⊥-2.414, .000⊥.000, 1.000⊥-.414, .000⊥.000, 1.000⊥.414, .000⊥.000, 1.000⊥2.414, 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,
BBC BASIC
<lang bbcbasic> @% = &60A
DIM Complex{r#, i#} DIM in{(7)} = Complex{}, out{(7)} = Complex{} DATA 1, 1, 1, 1, 0, 0, 0, 0 PRINT "Input (real, imag):" FOR I% = 0 TO 7 READ in{(I%)}.r# PRINT in{(I%)}.r# "," in{(I%)}.i# NEXT PROCfft(out{()}, in{()}, 0, 1, DIM(in{()},1)+1) PRINT "Output (real, imag):" FOR I% = 0 TO 7 PRINT out{(I%)}.r# "," out{(I%)}.i# NEXT END DEF PROCfft(b{()}, o{()}, B%, S%, N%) LOCAL I%, t{} : DIM t{} = Complex{} IF S% < N% THEN PROCfft(o{()}, b{()}, B%, S%*2, N%) PROCfft(o{()}, b{()}, B%+S%, S%*2, N%) FOR I% = 0 TO N%-1 STEP 2*S% t.r# = COS(-PI*I%/N%) t.i# = SIN(-PI*I%/N%) PROCcmul(t{}, o{(B%+I%+S%)}) b{(B%+I% DIV 2)}.r# = o{(B%+I%)}.r# + t.r# b{(B%+I% DIV 2)}.i# = o{(B%+I%)}.i# + t.i# b{(B%+(I%+N%) DIV 2)}.r# = o{(B%+I%)}.r# - t.r# b{(B%+(I%+N%) DIV 2)}.i# = o{(B%+I%)}.i# - t.i# NEXT ENDIF ENDPROC DEF PROCcmul(c{},d{}) LOCAL r#, i# r# = c.r#*d.r# - c.i#*d.i# i# = c.r#*d.i# + c.i#*d.r# c.r# = r# c.i# = i# ENDPROC
</lang> Output:
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.41421 0, 0 1, -0.414214 0, 0 1, 0.414214 0, 0 1, 2.41421
C
Inplace FFT with O(n) memory usage. Note: array size is assumed to be power of 2 and not checked by code; you can just pad it with 0 otherwise. Also, complex
is C99 standard.<lang C>#include <stdio.h>
- include <math.h>
- include <complex.h>
double PI; typedef double complex cplx;
void _fft(cplx buf[], cplx out[], int n, int step) { if (step < n) { _fft(out, buf, n, step * 2); _fft(out + step, buf + step, n, step * 2);
for (int i = 0; i < n; i += 2 * step) { cplx t = cexp(-I * PI * i / n) * out[i + step]; buf[i / 2] = out[i] + t; buf[(i + n)/2] = out[i] - t; } } }
void fft(cplx buf[], int n) { cplx out[n]; for (int i = 0; i < n; i++) out[i] = buf[i];
_fft(buf, out, n, 1); }
int main() { PI = atan2(1, 1) * 4; cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};
void show(char * s) { printf(s); for (int i = 0; i < 8; i++) if (!cimag(buf[i])) printf("%g ", creal(buf[i])); else printf("(%g, %g) ", creal(buf[i]), cimag(buf[i])); }
show("Data: "); fft(buf, 8); show("\nFFT : ");
return 0; }</lang>Output<lang>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)</lang>
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 = std::polar(1.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"; } return 0;
}</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*/ creal[] fft(/*in*/ creal[] x) /*pure nothrow*/ {
int N = x.length; if (N <= 1) return x; auto ev = stride(x, 2).array().fft(); auto od = stride(x[1 .. $], 2).array().fft(); auto l = map!(k => ev[k] + expi(-2*PI*k/N) * od[k])(iota(N / 2)); auto r = map!(k => ev[k] - expi(-2*PI*k/N) * od[k])(iota(N / 2)); return l.chain(r).array();
}
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]
But in D2 built-in complex numbers will be deprecated.
GAP
<lang gap># Here an implementation with no optimization (O(n^2)).
- In GAP, E(n) = exp(2*i*pi/n), a primitive root of the unity.
Fourier := function(a) local n, z; n := Size(a); z := E(n); return List([0 .. n - 1], k -> Sum([0 .. n - 1], j -> a[j + 1]*z^(-k*j))); end;
InverseFourier := function(a) local n, z; n := Size(a); z := E(n); return List([0 .. n - 1], k -> Sum([0 .. n - 1], j -> a[j + 1]*z^(k*j)))/n; end;
Fourier([1, 1, 1, 1, 0, 0, 0, 0]);
- [ 4, 1-E(8)-E(8)^2-E(8)^3, 0, 1-E(8)+E(8)^2-E(8)^3,
- 0, 1+E(8)-E(8)^2+E(8)^3, 0, 1+E(8)+E(8)^2+E(8)^3 ]
InverseFourier(last);
- [ 1, 1, 1, 1, 0, 0, 0, 0 ]</lang>
Go
<lang go>package main
import (
"fmt" "math" "math/cmplx"
)
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 := cmplx.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)
Haskell
<lang haskell>import Data.Complex
-- Cooley-Tukey fft [] = [] fft [x] = [x] fft xs = zipWith (+) ys ts ++ zipWith (-) ys ts
where n = length xs ys = fft evens zs = fft odds (evens, odds) = split xs split [] = ([], []) split [x] = ([x], []) split (x:y:xs) = (x:xt, y:yt) where (xt, yt) = split xs ts = zipWith (\z k -> exp' k n * z) zs [0..] exp' k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)
main = mapM_ print $ fft [1,1,1,1,0,0,0,0]</lang>
And the output:
4.0 :+ 0.0 1.0 :+ (-2.414213562373095) 0.0 :+ 0.0 1.0 :+ (-0.4142135623730949) 0.0 :+ 0.0 0.9999999999999999 :+ 0.4142135623730949 0.0 :+ 0.0 0.9999999999999997 :+ 2.414213562373095
J
Based on j:Essays/FFT, with some simplifications -- sacrificing accuracy, optimizations and convenience which are not relevant to the task requirements, 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.
Here is a representation of an example which appears in some of the other implementations, here: <lang J> Re=: {.@+.@fft
Im=: {:@+.@fft M=: 4#1 0 M
1 1 1 1 0 0 0 0
Re M
4 1 0 1 0 1 0 1
Im M
0 2.41421 0 0.414214 0 _0.414214 0 _2.41421</lang>
Note that Re and Im are not functions of 1 and 0 but are functions of the complete sequence.
Also note that floating point values are being calculated to the full 53 bits of mantissa precision of 64 bit floating point numbers, but only the first six digits are displayed, by default. Finally, note that J uses a different character for negative sign than for subtraction, to eliminate ambiguity (is this a single list of numbers or are lists being subtracted?).
Liberty BASIC
<lang lb>
P =8 S =int( log( P) /log( 2) +0.9999)
Pi =3.14159265 R1 =2^S
R =R1 -1 R2 =div( R1, 2) R4 =div( R1, 4) R3 =R4 +R2
Dim Re( R1), Im( R1), Co( R3)
for N =0 to P -1 read dummy: Re( N) =dummy read dummy: Im( N) =dummy next N
data 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0
S2 =div( S, 2) S1 =S -S2 P1 =2^S1 P2 =2^S2
dim V( P1 -1) V( 0) =0 DV =1 DP =P1
for J =1 to S1 HA =div( DP, 2) PT =P1 -HA for I =HA to PT step DP V( I) =V( I -HA) +DV next I DV =DV +DV DP =HA next J
K =2 *Pi /R1
for X =0 to R4 COX =cos( K *X) Co( X) =COX Co( R2 -X) =0 -COX Co( R2 +X) =0 -COX next X
print "FFT: bit reversal"
for I =0 to P1 -1 IP =I *P2 for J =0 to P2 -1 H =IP +J G =V( J) *P2 +V( I) if G >H then temp =Re( G): Re( G) =Re( H): Re( H) =temp if G >H then temp =Im( G): Im( G) =Im( H): Im( H) =temp next J next I
T =1
for stage =0 to S -1 print " Stage:- "; stage D =div( R2, T) for Z =0 to T -1 L =D *Z LS =L +R4 for I =0 to D -1 A =2 *I *T +Z B =A +T F1 =Re( A) F2 =Im( A) P1 =Co( L) *Re( B) P2 =Co( LS) *Im( B) P3 =Co( LS) *Re( B) P4 =Co( L) *Im( B) Re( A) =F1 +P1 -P2 Im( A) =F2 +P3 +P4 Re( B) =F1 -P1 +P2 Im( B) =F2 -P3 -P4 next I next Z T =T +T next stage
print " M Re( M) Im( M)"
for M =0 to R if abs( Re( M)) <10^-5 then Re( M) =0 if abs( Im( M)) <10^-5 then Im( M) =0 print " "; M, Re( M), Im( M) next M
end
wait
function div( a, b) div =int( a /b) end function
end
</lang>
M Re( M) Im( M) 0 4 0 1 1.0 -2.41421356 2 0 0 3 1.0 -0.41421356 4 0 0 5 1.0 0.41421356 6 0 0 7 1.0 2.41421356
Mathematica
Mathematica has a built-in FFT function which uses a proprietary algorithm developed at Wolfram Research. It also has an option to tune the algorithm for specific applications. The options shown below, while not default, produce output that is consistent with most other FFT routines.
<lang Mathematica> Fourier[{1,1,1,1,0,0,0,0}, FourierParameters->{1,-1}] </lang>
Output:
{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}
MATLAB / Octave
Matlab/Octave have a builtin FFT function.
<lang MATLAB> fft([1,1,1,1,0,0,0,0]') </lang> Output:
ans = 4.00000 + 0.00000i 1.00000 - 2.41421i 0.00000 + 0.00000i 1.00000 - 0.41421i 0.00000 + 0.00000i 1.00000 + 0.41421i 0.00000 - 0.00000i 1.00000 + 2.41421i
OCaml
This is a simple implementation of the Cooley-Tukey pseudo-code <lang OCaml>open Complex
let fac k n =
let m2pi = -4.0 *. acos 0.0 in polar 1.0 (m2pi*.(float k)/.(float n))
let merge l r n =
let f (k,t) x = (succ k, (mul (fac k n) x) :: t) in let z = List.rev (snd (List.fold_left f (0,[]) r)) in (List.map2 add l z) @ (List.map2 sub l z)
let fft lst =
let rec ditfft2 a n s = if n = 1 then [List.nth lst a] else let odd = ditfft2 a (n/2) (2*s) in let even = ditfft2 (a+s) (n/2) (2*s) in merge odd even n in ditfft2 0 (List.length lst) 1;;
let show l =
let pr x = Printf.printf "(%f %f) " x.re x.im in (List.iter pr l; print_newline ()) in
let indata = [one;one;one;one;zero;zero;zero;zero] in show indata; show (fft indata)</lang>
Output:
(1.000000 0.000000) (1.000000 0.000000) (1.000000 0.000000) (1.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (0.000000 0.000000) (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)
PARI/GP
Naive implementation, using the same testcase as Ada: <lang parigp>FFT(v)=my(t=-2*Pi*I/#v,tt);vector(#v,k,tt=t*(k-1);sum(n=0,#v-1,v[n+1]*exp(tt*n))); FFT([1,1,1,1,0,0,0,0])</lang> Output:
[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]
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
PicoLisp
<lang PicoLisp># apt-get install libfftw3-dev
(scl 4)
(de FFTW_FORWARD . -1) (de FFTW_ESTIMATE . 64)
(de fft (Lst)
(let (Len (length Lst) In (native "libfftw3.so" "fftw_malloc" 'N (* Len 16)) Out (native "libfftw3.so" "fftw_malloc" 'N (* Len 16)) P (native "libfftw3.so" "fftw_plan_dft_1d" 'N Len In Out FFTW_FORWARD FFTW_ESTIMATE ) ) (struct In NIL (cons 1.0 (apply append Lst))) (native "libfftw3.so" "fftw_execute" NIL P) (prog1 (struct Out (make (do Len (link (1.0 . 2))))) (native "libfftw3.so" "fftw_destroy_plan" NIL P) (native "libfftw3.so" "fftw_free" NIL Out) (native "libfftw3.so" "fftw_free" NIL In) ) ) )</lang>
Test: <lang PicoLisp>(for R (fft '((1.0 0) (1.0 0) (1.0 0) (1.0 0) (0 0) (0 0) (0 0) (0 0)))
(tab (6 8) (round (car R)) (round (cadr R)) ) )</lang>
Output:
4.000 0.000 1.000 -2.414 0.000 0.000 1.000 -0.414 0.000 0.000 1.000 0.414 0.000 0.000 1.000 2.414
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
REXX
This REXX program is modeled after the Run BASIC version and is a radix-2 DIC (decimation-in-time) form of the
Cooley-Turkey FFT algorithm, and as such, this simplified form assumes that the number of data points is equal to an exact power of two.
Note that the REXX language doesn't have any higher math functions, so the functions COS and R2R
(normalization of radians) have been included here, as well as the constant pi .
This program also adds zero values if the number of data points in the list doesn't exactly equal to a power of two.
This is called zero-padding.
<lang rexx>/*REXX pgm does a fast Fourier transform (FFT) on a set of complex nums.*/
numeric digits 80
arg data; if data= then data='1 1 1 1 0' /*no data? Then use default*/
data=translate(data,'J',"I") /*allow use of i as well as j */
size=words(data)
do sig=0 until 2**sig>=size; end /*# of args exactly a power of 2?*/ do j=size+1 to 2**sig;data=data 0;end /*add zeroes until a power of 2*/
size=words(data); call hd /*┌─────────────────────────────┐*/
/*│ Numbers in data can be in │*/ do j=0 for size /*│ 7 formats: real │*/ _=word(data,j+1) /*│ real,imag │*/ parse var _ #.1.j ',' #.2.j /*│ ,imag │*/ if right(#.1.j,1)=='J' then /*│ nnnJ │*/ parse var #.1.j #2.j "J" @.1.j /*│ nnnj │*/ do p=1 for 2 /*omitted?*/ /*│ nnnI │*/ #.p.j=word(#.p.j 0,1) /*│ nnni │*/ end /*└─────────────────────────────┘*/ say "FFT in " center(j+1,3) ' ' nice(#.1.j) nice(#.2.j,'i') end
say; tran=2*pi()/2**sig; !.=0 hsig=2**sig%2; counterA=2**(sig-sig%2); pointer=counterA; doubler=1
do sig-sig%2; halfpointer=pointer%2 do i=halfpointer by pointer to counterA-halfpointer _=i-halfpointer; !.i=!._+doubler end doubler=doubler*2; pointer=halfpointer end do j=0 to 2**sig%4; cmp.j=cos(j*tran); _m=hsig-j; _p=hsig+j cmp._m=-cmp.j; cmp._p=-cmp.j end
counterB=2**(sig%2)
do i=0 for counterA; p=i*counterB do j=0 for counterB; h=p+j; _=!.j*counterB+!.i; if _<=h then iterate parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._ end /*(above) switch 2 sets of values*/ end
double=1
do sig; w=2**sig%2%double do k=0 for double; lb=w*k; lh=lb+2**sig%4 do j=0 for w; a=j*double*2+k; b=a+double; r=#.1.a; i=#.2.a c1=cmp.lb*#.1.b; c4=cmp.lb*#.2.b c2=cmp.lh*#.2.b; c3=cmp.lh*#.1.b #.1.a=r+c1-c2; #.2.a=i+c3+c4 #.1.b=r-c1+c2; #.2.b=i-c3-c4 end end double=double*2 end
call hd; do i=0 for size
say "FFT out" center(i+1,3) ' ' nice(#.1.i) nice(#.2.i,'j') end
exit /*─────────────────────────────────────subroutines──────────────────────*/ hd:_='numbers # real-part imaginary-part';say _;say translate(_," "copies('═',255),xrange());return pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862;return pi r2r: return arg(1)//(2*pi()) /*reduce radians to unit circle. */ cos: procedure; arg x; x=r2r(x); return .sincos(1,1,-1) .sincos: parse arg z,_,i; x=x*x; p=z
do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end; return z /*subroutine to make complex numbers look nicer for output.*/
nice: procedure; parse arg x,j; numeric digits digits()%10; nz='1e-'digits()
if abs(x)<nz then x=0;x=x/1;if x=0 & j\== then return x=format(x,,digits()); if pos('.',x)\==0 then x=strip(x,'T',0) x=strip(x,,'.'); if x>=0 then x=' '||x;return left(x||j,digits()+4)</lang>
Output when using the default input of: 1 1 1 1 0
numbers # real-part imaginary-part ═══════ ═ ═════════ ══════════════ FFT in 1 1 FFT in 2 1 FFT in 3 1 FFT in 4 1 FFT in 5 0 FFT in 6 0 FFT in 7 0 FFT in 8 0 numbers # real-part imaginary-part ═══════ ═ ═════════ ══════════════ FFT out 1 4 FFT out 2 1 -2.4142136j FFT out 3 0 FFT out 4 1 -0.41421356j FFT out 5 0 FFT out 6 1 0.41421356j FFT out 7 0 FFT out 8 1 2.4142136j
Run BASIC
<lang runbasic>cnt = 8 sig = int(log(cnt) /log(2) +0.9999)
pi = 3.14159265 real1 = 2^sig
real = real1 -1 real2 = int(real1 / 2) real4 = int(real1 / 4) real3 = real4 +real2
dim rel(real1) dim img(real1) dim cmp(real3)
for i = 0 to cnt -1
read rel(i) read img(i)
next i
data 1,0, 1,0, 1,0, 1,0, 0,0, 0,0, 0,0, 0,0
sig2 = int(sig / 2) sig1 = sig -sig2 cnt1 = 2^sig1 cnt2 = 2^sig2
dim v(cnt1 -1) v(0) = 0 dv = 1 ptr = cnt1
for j = 1 to sig1
hlfPtr = int(ptr / 2) pt = cnt1 - hlfPtr for i = hlfPtr to pt step ptr v(i) = v(i -hlfPtr) + dv next i dv = dv + dv ptr = hlfPtr
next j
k = 2 *pi /real1
for x = 0 to real4
cmp(x) = cos(k *x) cmp(real2 - x) = 0 - cmp(x) cmp(real2 + x) = 0 - cmp(x)
next x
print "fft: bit reversal"
for i = 0 to cnt1 -1
ip = i *cnt2 for j = 0 to cnt2 -1 h = ip +j g = v(j) *cnt2 +v(i) if g >h then temp = rel(g) rel(g) = rel(h) rel(h) = temp temp = img(g) img(g) = img(h) img(h) = temp end if next j
next i
t = 1 for stage = 1 to sig
print " stage:- "; stage d = int(real2 / t) for ii = 0 to t -1 l = d *ii ls = l +real4 for i = 0 to d -1 a = 2 *i *t +ii b = a +t f1 = rel(a) f2 = img(a) cnt1 = cmp(l) *rel(b) cnt2 = cmp(ls) *img(b) cnt3 = cmp(ls) *rel(b) cnt4 = cmp(l) *img(b) rel(a) = f1 + cnt1 - cnt2 img(a) = f2 + cnt3 + cnt4 rel(b) = f1 - cnt1 + cnt2 img(b) = f2 - cnt3 - cnt4 next i next ii t = t +t
next stage
print " Num real imag" for i = 0 to real
if abs(rel(i)) <10^-5 then rel(i) = 0 if abs(img(i)) <10^-5 then img(i) = 0 print " "; i;" ";using("##.#",rel(i));" ";img(i)
next i end</lang>
Num real imag 0 4.0 0 1 1.0 -2.41421356 2 0.0 0 3 1.0 -0.414213565 4 0.0 0 5 1.0 0.414213562 6 0.0 0 7 1.0 2.41421356
Scala
</lang>
Tcl
<lang tcl>package require math::constants package require math::fourier
math::constants::constants pi
- 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
- Uses FFT if input length is power of 2, and a less efficient algorithm otherwise
set fft [math::fourier::dft $wave]
- 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
Ursala
The fftw
library is callable from Ursala using the syntax ..u_fw_dft
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.
<lang ursala>#import nat
- import flo
f = <1+0j,1+0j,1+0j,1+0j,0+0j,0+0j,0+0j,0+0j> # complex sequence of 4 1's and 4 0's
g = c..mul^*D(sqrt+ float+ length,..u_fw_dft) f # its fft
- cast %jLW
t = (f,g)</lang> output:
( < 1.000e+00+0.000e+00j, 1.000e+00+0.000e+00j, 1.000e+00+0.000e+00j, 1.000e+00+0.000e+00j, 0.000e+00+0.000e+00j, 0.000e+00+0.000e+00j, 0.000e+00+0.000e+00j, 0.000e+00+0.000e+00j>, < 4.000e+00+0.000e+00j, 1.000e+00-2.414e+00j, 0.000e+00+0.000e+00j, 1.000e+00-4.142e-01j, 0.000e+00+0.000e+00j, 1.000e+00+4.142e-01j, 0.000e+00+0.000e+00j, 1.000e+00+2.414e+00j>)