Anonymous user
Fast Fourier transform: Difference between revisions
m
→{{header|REXX}}: added/changed comments and whitespace, changed indentations.
(→{{header|C}}: added apple platform library usage) |
m (→{{header|REXX}}: added/changed comments and whitespace, changed indentations.) |
||
Line 2,202:
<br>('''cos'''ine and reduce radians to a unit circle).
A normalization of radians function  ' ('''r2r''') has been included here, as well as the constant '''pi'''.
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''.
<lang rexx>/*REXX
numeric digits length( pi() ) - 1 /*limited by the PI function result. */
arg data /*ARG verb uppercases the DATA from CL.*/
if data='' then data=1 1 1 1 0 /*Not specified? Then use the default.*/
size=words(data); pad=left('',6) /*PAD: for indenting and padding SAYs.*/
do p=0 until 2**p>=size ; end /*number of args exactly a power of 2? */
do j=size+1 to 2**p;data=data 0; end /*add zeroes to DATA 'til a power of 2.*/
size=words(data); ph=p%2; call hdr /*
/* [↓] TRANSLATE allows I & J*/
do j=0 for size /*║ seven formats:
_=translate( word(data,j+1), 'J', "I") /*║
parse var _ #.1.j '' $ 1
if $=='J' then parse var #.1.j #2.j
do m=1 for 2; #.m.j= word(#.m.j 0, 1) /*║
end /*m*/ /*omitted part? [↑]
/*
say pad
end /*j*/
say
tran=pi()*2 / 2**p; !.=0;
say
do p-ph; halfPtr=ptr%2
do i=halfPtr by ptr to A-halfPtr; _=i-halfPtr; !.i=!._+dbl
end /*i*/
dbl=dbl*2; ptr=halfPtr
end /*p-ph*/
B=2**ph
end /*j*/ /* [↑] swap two sets of values. */
dbl=1
do p ; w=hp % dbl
r=#.1.a; i=#.2.a ;
#.1.a=r + c1 - c2
end
end
end
▲ end /*p*/
call hdr
do i=0 for size
say pad " FFT out " center(i+1,7) pad fmt(#.1.i) fmt(#.2.i,'j')
end /*i*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
cos: procedure; parse arg x; q=r2r(x)**2; z=1; _=1;
do k=2 by 2; _=-_*q/(k*(k-1)); z=z+_; if z=p then leave; p=z; end; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
fmt: procedure; parse arg y,j; y=y/1 /*
if abs(y) < '1e-'digits()%4 then y=0; if y=0 & j\=='' then return ''
y=format(y, , 10); if pos(.,y)\==0 then y=strip(y, 'T', 0)
y=strip(y, , .);
/*──────────────────────────────────────────────────────────────────────────────────────*/
hdr: _='───data─── num real─part imaginary─part'; say pad _
say pad translate(_, " "copies('═',256), " "xrange()); return
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
/*──────────────────────────────────────────────────────────────────────────────────────*/
r2r: return arg(1) // ( pi()*2 ) /*reduce the radians to a unit circle. */</lang>
'''output''' when using the default input of: <tt> 1 1 1 1 0 </tt>
<pre>
| |||