Anonymous user
Fast Fourier transform: Difference between revisions
→{{header|REXX}}: added the REXX language. -- ~~~~
(Added PicoLisp) |
(→{{header|REXX}}: added the REXX language. -- ~~~~) |
||
Line 761:
Output:
<pre>4+0.000000i 1-2.414214i 0+0.000000i 1-0.414214i 0+0.000000i 1+0.414214i 0+0.000000i 1+2.414214i</pre>
=={{header|REXX}}==
This REXX program is modeled after the <tt> Run BASIC </tt>
<br><br>Note that the REXX language doesn't have any higher math functions, so the functions
<br>LN (natural log), COS, and normalization of radians have been included here,
<br>as well as two constants (e and pi).
<lang rexx>
/*REXX pgm does a fast Fourier transform (FFT) on a set of complex nums.*/
numeric digits 80 /*unknown input, so use high prec*/
arg data /*┌─────────────────────────────┐*/
if data='' then data=1 1 1 1 0 0 0 0 /*│ Numbers in data can be in │*/
size=words(data) /*│ two formats: │*/
call hdr /*│ real │*/
do j=0 for size /*│ real,imsg │*/
_=word(data,j+1) /*└─────────────────────────────┘*/
parse var _ #.1.j ',' #.2.j
do p=1 for 2; if #.p.j=='' then #.p.j=0; end /*handle omits*/
say "FFT in " center(j+1,3) ' ' nice(#.1.j) nice(#.2.j,'i')
end
sig=1+trunc(ln(size)/ln(2)); call pi; tran=2*pi/2**sig; !.=0
hsig=2**sig%2; counterA=2**(sig-sig%2); pointer=counterA; doubler=1
do sig-sig%2; halfpointer=pointer%2
do i=halfpointer by pointer to counterA-halfpointer
_=i-halfpointer; !.i=!._+doubler
end
doubler=doubler*2; pointer=halfpointer
end
do j=0 to 2**sig%4; cmp.j=cos(j*tran); _m=hsig-j; _p=hsig+j
cmp._m=-cmp.j; cmp._p=-cmp.j
end
counterB=2**(sig%2)
do i=0 for counterA; p=i*counterB
do j=0 for counterB; h=p+j; _=!.j*counterB+!.i; if _<=h then iterate
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._
end /*(above) switch 2 sets of values*/
end
/*─────────────────────────────────────perform FFT transformation of #s.*/
double=1
do sig; w=2**sig%2%double
do k=0 for double; lb=w*k; lh=lb+2**sig%4
do j=0 for w; a=j*double*2+k; b=a+double; 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
#.1.a=r+c1-c2; #.2.a=i+c3+c4
#.1.b=r-c1+c2; #.2.b=i-c3-c4
end
end
double=double*2
end
call hdr; do i=0 for size
say "FFT out" center(i+1,3) ' ' nice(#.1.i) nice(#.2.i,'i')
end
exit
/*─────────────────────────────────────subroutines──────────────────────*/
hdr: say; say 'numbers # real-part imaginary-part'
say '─────── ─── ───────── ──────────────' ;return
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862;return pi
e: e=2.7182818284590452353602874713526624977572470936999595749669676277240766303535;return e
r2r: return arg(1)//(2*pi()) /*reduce radians to unit circle. */
cos: procedure; arg x; x=r2r(x); return .sincos(1,1,-1)
.sincos: parse arg z,_,i; x=x*x; p=z
do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end; return z
/*subroutine to make complex numbers look nicer for output.*/
nice: procedure; parse arg x,j; numeric digits digits()%10; nz='1e-'digits()
if abs(x)<nz then x=0;x=x/1;if x=0 & j\=='' then return ''
x=format(x,,digits()); if pos('.',x)\==0 then x=strip(x,'T',0)
x=strip(x,,'.'); if x>=0 then x=' '||x;return left(x||j,digits()+4)
ln: procedure; parse arg x,f; call e; ig=x>1.5; is=1-2*(ig\==1); ii=0; xx=x; return .ln_comp()
.ln_comp: do while ig & xx>1.5 | \ig & xx<.5; _=e; do k=-1; iz=xx*_**-is;
if k>=0&(ig&iz<1|\ig&iz>.5) then leave; _=_*_;izz=iz;end; xx=izz; ii=ii+is*2**k; end
x=x*e**-ii-1; z=0; _=-1; p=z; do k=1; _=-_*x; z=z+_/k; if z=p then leave; p=z; end
return z+ii
</lang>
Output when using the default input:
<pre style="height:60ex;overflow:scroll">
numbers # real-part imaginary-part
─────── ─── ───────── ──────────────
FFT in 1 1
FFT in 2 1
FFT in 3 1
FFT in 4 1
FFT in 5 0
FFT in 6 0
FFT in 7 0
FFT in 8 0
numbers # real-part imaginary-part
─────── ─── ───────── ──────────────
FFT out 1 4
FFT out 2 1 -2.4142136i
FFT out 3 0
FFT out 4 1 -0.41421356i
FFT out 5 0
FFT out 6 1 0.41421356i
FFT out 7 0
FFT out 8 1 2.4142136i
</pre>
=={{header|Run BASIC}}==
| |||