Fast Fourier transform: Difference between revisions

Content added Content deleted
m (→‎{{header|ALGOL 68}}: misc code refactor)
(→‎{{header|ALGOL 68}}: replace array SLICE operator with array DICE operator... simpler.)
Line 97: Line 97:
=={{header|ALGOL 68}}==
=={{header|ALGOL 68}}==
{{trans|Python}}Note: This specimen retains the original [[#Python|Python]] coding style.
{{trans|Python}}Note: This specimen retains the original [[#Python|Python]] coding style.
{{works with|ALGOL 68|Revision 1 - no extension to language.}}
{{works with|ALGOL 68|Revision 1 - one minor extension to language used - PRAGMA READ, similar to C's #include directive.}}
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-2.3.5 algol68g-2.3.5].}}
{{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''.}}
{{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: Fast_Fourier_transform.a68'''<lang algol68>#!/usr/local/bin/a68g --script #
'''File: Template.Fast_Fourier_transform.a68'''<lang algol68>PRIO DICE = 9; # ideally = 11 #
⚫
# -*- coding: utf-8 -*- #

⚫
MODE VALUE = COMPL;

MODE OPTINT = UNION(INT, VOID),
SLICE = STRUCT(OPTINT lwb, upb, step);

⚫
FORMAT slice fmt = $"["g(0)":"g(0)"/"g(0)"] "$;

OP ELEM = (SLICE slice, []VALUE in)[]VALUE: (
INT lwb = (lwb OF slice|(INT lwb ): lwb |LWB in);
INT upb = (upb OF slice|(INT upb ): upb |UPB in);
INT step = (step OF slice|(INT step): step|1);


OP DICE = ([]SCALAR in, INT step)[]SCALAR: (
### Dice the array, extract array values a "step" apart ###
IF step = 1 THEN
IF step = 1 THEN
in[lwb:upb]
in
ELSE
ELSE
⚫
[(upb-lwb)%step+1]VALUE out;
INT upb out := 0;
INT upb out := 0;
⚫
[(UPB in-LWB in)%step+1]SCALAR out;

FOR index FROM lwb BY step TO upb DO
FOR index FROM LWB in BY step TO UPB in DO
out[upb out+:=1] := in[index] OD;
out[upb out+:=1] := in[index] OD;
out
out[@LWB in]
FI
FI
);
);


PROC fft = ([]VALUE in)[]VALUE: (
PROC fft = ([]SCALAR in t)[]SCALAR: (
# ASSERT LWB in = 1 required #
### The Cooley-Tukey FFT algorithm ###
IF LWB in >= UPB in THEN
IF LWB in t >= UPB in t THEN
in
in t[@0]
ELSE
ELSE
INT n = UPB in - LWB in + 1;
[]SCALAR t = in t[@0];
[LWB in: UPB in]VALUE out;
INT n = UPB t + 1, half n = n % 2;
[LWB t:UPB t]SCALAR coef;
[]VALUE even = fft(SLICE(EMPTY, EMPTY,2) ELEM in)[@LWB in],
odd = fft(SLICE(LWB in+1,EMPTY,2) ELEM in)[@LWB in];
INT mid in = (UPB in + LWB in) % 2;


FOR i FROM LWB in TO mid in DO
[]SCALAR even = fft(t DICE 2),
INT k = i - LWB in;
odd = fft(t[1:]DICE 2);

⚫
COMPL exp = complex exp(-(0I2)*pi*k/n)*odd[i];
out[i] := even[i] + exp;
COMPL i = 0 I 1;

⚫
out[mid in + i] := even[i] - exp
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;
OD;
out
coef
FI
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,real width % 3-2)$;
FORMAT real fmt := $g(0,3)$;
FORMAT compl fmt = $f(real fmt)"⊥"f(real fmt)$;
FORMAT real array fmt := $f(real fmt)", "$;
FORMAT array compl fmt = $f(compl fmt)", "$;
FORMAT compl fmt := $f(real fmt)"⊥"f(real fmt)$;
⚫
FORMAT compl array fmt := $f(compl fmt)", "$;


test:(
⚫
printf((
[]COMPL
$f(array compl fmt)$, fft((1, 1, 1, 1, 0, 0, 0, 0)), $l$,
tooth wave ft = fft((1, 1, 1, 1, 0, 0, 0, 0)),
$f(array compl fmt)$,
fft((0, 0.924, 0.707,-0.383,-1,-0.383, 0.707, 0.924,
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)), $l$
0,-0.924,-0.707, 0.383, 1, 0.383,-0.707,-0.924));
⚫
printf((
⚫
))</lang>'''Output:'''
$"Tooth wave: "$,compl array fmt, tooth wave ft, $l$,
$"1¼ cycle wave: "$, compl array fmt, one and a quarter wave ft, $l$
))
⚫
)</lang>'''Output:'''
<pre>
<pre>
4.000⊥.000, 1.000⊥-2.414, .000⊥.000, 1.000⊥-.414, .000⊥.000, 1.000⊥.414, .000⊥.000, 1.000⊥2.414,
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,
.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,
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>
</pre>


Retrieved from "https://rosettacode.org/wiki/Fast_Fourier_transform"