Fast Fourier transform: Difference between revisions
Content added Content deleted
(→{{header|Perl}}: fix missing actual function call) |
m (→{{header|REXX}}: changed whitespace and indentations, optimized a computation.) |
||
| Line 2,382: | Line 2,382: | ||
arg data /*ARG verb uppercases the DATA from CL.*/ |
arg data /*ARG verb uppercases the DATA from CL.*/ |
||
if data='' then data=1 1 1 1 0 /*Not specified? Then use the default.*/ |
if data='' then data=1 1 1 1 0 /*Not specified? Then use the default.*/ |
||
size=words(data); pad=left('',6) |
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 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.*/ |
do j=size+1 to 2**p; data=data 0; end /*add zeroes to DATA 'til a power of 2.*/ |
||
| Line 2,398: | Line 2,398: | ||
end /*j*/ |
end /*j*/ |
||
say |
say |
||
tran=pi()*2 / 2**p; !.=0; hp=2**p %2; A=2**(p-ph); ptr=A; dbl=1 |
tran=pi()*2 / 2**p; !.=0; hp=2**p %2; A=2**(p-ph); ptr=A; dbl=1 |
||
say |
say |
||
do p-ph; halfPtr=ptr % 2 |
do p-ph; halfPtr=ptr % 2 |
||
do i=halfPtr by ptr to A-halfPtr; |
do i=halfPtr by ptr to A-halfPtr; _=i - halfPtr; !.i=!._ + dbl |
||
end /*i*/ |
end /*i*/ |
||
ptr=halfPtr; dbl=dbl+dbl |
|||
end /*p-ph*/ |
end /*p-ph*/ |
||
do j=0 to 2**p%4; cmp.j=cos(j*tran); _=hp - j; cmp._= -cmp.j |
do j=0 to 2**p%4; cmp.j=cos(j*tran); _=hp - j; cmp._= -cmp.j |
||
_=hp + j; cmp._= -cmp.j |
_=hp + j; cmp._= -cmp.j |
||
end /*j*/ |
end /*j*/ |
||
B=2**ph |
B=2**ph |
||
do i=0 for A; q=i * B |
do i=0 for A; q=i * B |
||
do j=0 for B; h=q+j; _=!.j*B+!.i; if _<=h then iterate |
do j=0 for B; h=q+j; _=!.j*B+!.i; if _<=h then iterate |
||
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._ |
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._ |
||
end /*j*/ /* [↑] swap two sets of values. */ |
end /*j*/ /* [↑] swap two sets of values. */ |
||
end /*i*/ |
end /*i*/ |
||
dbl=1 |
dbl=1 |
||
do p ; w=hp % dbl |
do p ; w= hp % dbl |
||
do k=0 for dbl ; Lb=w * k ; Lh=Lb + 2**p % 4 |
do k=0 for dbl ; Lb= w * k ; Lh= Lb + 2**p % 4 |
||
do j=0 for w ; a=j * dbl * 2 + k ; b= a + dbl |
do j=0 for w ; a= j * dbl * 2 + k ; b= a + dbl |
||
r=#.1.a; i=#.2.a ; c1=cmp.Lb * #.1.b ; c4=cmp.Lb * #.2.b |
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 |
c2= cmp.Lh * #.2.b ; c3= cmp.Lh * #.1.b |
||
#.1.a=r + c1 - c2 ; #.2.a=i + c3 + c4 |
#.1.a= r + c1 - c2 ; #.2.a= i + c3 + c4 |
||
#.1.b=r - c1 + c2 ; #.2.b=i - c3 - c4 |
#.1.b= r - c1 + c2 ; #.2.b= i - c3 - c4 |
||
end /*j*/ |
end /*j*/ |
||
end /*k*/ |
end /*k*/ |
||
| Line 2,429: | Line 2,429: | ||
call hdr |
call hdr |
||
do i=0 for size |
do i=0 for size |
||
say pad " FFT out " center(i+1,7) pad fmt(#.1.i) fmt(#.2.i,'j') |
say pad " FFT out " center(i+1,7) pad fmt(#.1.i) fmt(#.2.i,'j') |
||
end /*i*/ /*[↑] #s are shown with 10 decimal digs*/ |
end /*i*/ /*[↑] #s are shown with 10 decimal digs*/ |
||
exit /*stick a fork in it, we're all done. */ |
exit /*stick a fork in it, we're all done. */ |
||
| Line 2,437: | Line 2,437: | ||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
fmt: procedure; parse arg y,j; y=y/1 /*prettifies complex numbers for output*/ |
fmt: procedure; parse arg y,j; y=y/1 /*prettifies complex numbers for output*/ |
||
if abs(y) < '1e-'digits()%4 then y=0; if y=0 & j\=='' then return '' |
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=format(y, , 10); if pos(.,y)\==0 then y=strip(y, 'T', 0) |
||
y=strip(y, , .); if y>=0 then y=' 'y; |
y=strip(y, , .); if y>=0 then y=' 'y; return left(y || j, 12) |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
hdr: _= '───data─── num real─part imaginary─part'; say pad _ |
hdr: _= '───data─── num real─part imaginary─part'; say pad _ |
||
| Line 2,446: | Line 2,446: | ||
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862 |
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862 |
||
r2r: return arg(1) // ( pi()*2 ) /*reduce the radians to a unit circle. */</lang> |
r2r: return arg(1) // ( pi()*2 ) /*reduce the radians to a unit circle. */</lang> |
||
{{out|output|text= when using the default input of: <tt> 1 1 1 1 0 </tt>}} |
|||
<pre> |
<pre> |
||
───data─── num real─part imaginary─part |
───data─── num real─part imaginary─part |
||