Fast Fourier transform: Difference between revisions

→‎{{header|ALGOL 68}}: replace array SLICE operator with array DICE operator... simpler.
m (→‎{{header|ALGOL 68}}: misc code refactor)
(→‎{{header|ALGOL 68}}: replace array SLICE operator with array DICE operator... simpler.)
Line 97:
=={{header|ALGOL 68}}==
{{trans|Python}}Note: This specimen retains the original [[#Python|Python]] coding style.
{{works with|ALGOL 68|Revision 1 - noone 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].}}
{{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'''<lang algol68>PRIO DICE = 9; #!/usr/local/bin/a68g --scriptideally = 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
in[lwb:upb]
ELSE
[(upb-lwb)%step+1]VALUE out;
INT upb out := 0;
[(upbUPB in-lwbLWB in)%step+1]VALUESCALAR out;
 
FOR index FROM lwbLWB in BY step TO upbUPB in DO
out[upb out+:=1] := in[index] OD;
out[@LWB in]
FI
);
 
PROC fft = ([]VALUESCALAR in t)[]VALUESCALAR: (
### ASSERTThe LWBCooley-Tukey inFFT = 1 requiredalgorithm ###
IF LWB in t >= UPB in t THEN
in t[@0]
ELSE
INT[]SCALAR nt = UPB in - LWB in + 1t[@0];
[LWBINT in:n = UPB in]VALUEt + 1, half n = n % out2;
[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[]SCALAR ieven FROM= LWBfft(t in TO mid inDICE DO2),
INT k = i - LWB inodd = fft(t[1:]DICE 2);
 
COMPL exp = complex exp(-(0I2)*pi*k/n)*odd[i];
COMPL out[i] := even[i]0 +I exp1;
 
out[mid in + i] := even[i] - exp
REAL w = 2*pi / n;
FOR k FROM LWB t TO half n-1 DO
COMPL expcis t = complexscalar exp(-0 I (0I2)*pi-w * k/n))*odd[ik];
coef[k] := even[k] + cis t;
outcoef[mid ink + ihalf n] := even[ik] - expcis t
OD;
outcoef
FI
);</lang>'''File: test.Fast_Fourier_transform.a68'''<lang algol68>#!/usr/local/bin/a68g --script #
);
# -*- coding: utf-8 -*- #
 
MODE VALUESCALAR = 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 complreal array fmt := $f(real fmt)", "f(real fmt)$;
FORMAT array compl fmt := $f(complreal fmt)", "f(real fmt)$;
FORMAT slicecompl array fmt := $"["gf(0)":"g(0)"/"g(0compl fmt)"], "$;
 
test:(
printf((
[]COMPL
$f(array compl fmt)$,tooth wave ft = fft((1, 1, 1, 1, 0, 0, 0, 0)), $l$,
$f(array compl fmt)$,
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$;
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>
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,
</pre>