Anonymous user
Fast Fourier transform: Difference between revisions
m
→{{header|REXX}}: added/changed comments and whitespace, showed more decimal digits in the output, simplified the FMT function.
(Adding PHP FFT & iFFT example) |
m (→{{header|REXX}}: added/changed comments and whitespace, showed more decimal digits in the output, simplified the FMT function.) |
||
Line 2,772:
<br>power of two. This is known as ''zero-padding''.
<lang rexx>/*REXX program performs a fast Fourier transform (FFT) on a set of complex numbers. */
numeric digits length( pi() ) -
arg data /*ARG verb uppercases the DATA from CL.*/
if data='' then data= 1 1 1 1 0
size=words(data); pad= left('',
do p=0 until 2**p>=size ; end
do j=size+1 to 2**p; data= data 0; end
size= words(data);
/* [↓] TRANSLATE allows I & J*/ /*║ Numbers in data can be in ║*/
do j=0 for size /*║ seven formats: real ║*/
_= translate( word(data, j+1), 'J', "I")
parse var _ #.1.j '' $ 1 "," #.2.j /*║ ,imag ║*/
if $=='J' then parse var #.1.j #2.j "J" #.1.j /*║ nnnJ ║*/
Line 2,791:
end /*j*/
say
tran= pi()*2 / 2**p;
say
do p-ph; halfPtr=ptr % 2
ptr= halfPtr; dbl= dbl + dbl
end /*p-ph*/
do j=0 to 2**p%4; cmp.j= cos(j*tran); _= hp - j; cmp._= -cmp.j
_= hp + j; cmp._= -cmp.j
end /*j*/
B= 2**ph
do i=0 for A; q= i * B
do j=0 for B; h=q+j; _= !.j*B+!.i;
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._
end /*j*/ /* [↑] swap two sets of values. */
end /*i*/
dbl= 1
do p
do k=0 for dbl
do j=0 for w
r= #.1.a; i= #.2.a ;
c2= cmp.Lh * #.2.b ; c3= cmp.Lh * #.1.b
#.1.a= r + c1 - c2 ; #.2.a= i + c3 + c4
Line 2,818:
end /*j*/
end /*k*/
dbl= dbl + dbl
end /*p*/
call hdr
do
say pad " FFT out " center(
end /*
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
cos: procedure;
do k=2 by 2; _=-_*q/(k*(k-1)); z=z+_; if z=p then return z; p=z; end /*k*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
fmt: procedure; parse arg y,j;
if abs(y) < '1e-'digits() %4 then y= 0;
dp= digits()%3; y= format(y, dp%6+1,
y= strip(y, 'T', .);
/*──────────────────────────────────────────────────────────────────────────────────────*/
hdr: _=pad '
say
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
r2r: return arg(1) // ( pi() * 2 )
{{out|output|text= when using the default
<pre>
</pre>
| |||