Another update from ingydotnet^djgoku
This commit is contained in:
Ingy döt Net
parent
91df62d461
commit
948b86eafa
7604 changed files with 108452 additions and 22726 deletions
|
|
@ -7,7 +7,8 @@ const double PI = 3.141592653589793238460;
|
|||
typedef std::complex<double> Complex;
|
||||
typedef std::valarray<Complex> CArray;
|
||||
|
||||
// Cooley–Tukey FFT (in-place)
|
||||
// Cooley–Tukey FFT (in-place, divide-and-conquer)
|
||||
// Higher memory requirements and redundancy although more intuitive
|
||||
void fft(CArray& x)
|
||||
{
|
||||
const size_t N = x.size();
|
||||
|
|
@ -30,6 +31,56 @@ void fft(CArray& x)
|
|||
}
|
||||
}
|
||||
|
||||
// Cooley-Tukey FFT (in-place, breadth-first, decimation-in-frequency)
|
||||
// Better optimized but less intuitive
|
||||
void fft(CArray &x)
|
||||
{
|
||||
// DFT
|
||||
unsigned int N = x.size(), k = N, n;
|
||||
double thetaT = 3.14159265358979323846264338328L / N;
|
||||
Complex phiT = Complex(cos(thetaT), sin(thetaT)), T;
|
||||
while (k > 1)
|
||||
{
|
||||
n = k;
|
||||
k >>= 1;
|
||||
phiT = phiT * phiT;
|
||||
T = 1.0L;
|
||||
for (unsigned int l = 0; l < k; l++)
|
||||
{
|
||||
for (unsigned int a = l; a < N; a += n)
|
||||
{
|
||||
unsigned int b = a + k;
|
||||
Complex t = x[a] - x[b];
|
||||
x[a] += x[b];
|
||||
x[b] = t * T;
|
||||
}
|
||||
T *= phiT;
|
||||
}
|
||||
}
|
||||
// Decimate
|
||||
unsigned int m = (unsigned int)log2(N);
|
||||
for (unsigned int a = 0; a < N; a++)
|
||||
{
|
||||
unsigned int b = a;
|
||||
// Reverse bits
|
||||
b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
|
||||
b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
|
||||
b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
|
||||
b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
|
||||
b = ((b >> 16) | (b << 16)) >> (32 - m);
|
||||
if (b > a)
|
||||
{
|
||||
Complex t = x[a];
|
||||
x[a] = x[b];
|
||||
x[b] = t;
|
||||
}
|
||||
}
|
||||
//// Normalize (This section make it not working correctly)
|
||||
//Complex f = 1.0 / sqrt(N);
|
||||
//for (unsigned int i = 0; i < N; i++)
|
||||
// x[i] *= f;
|
||||
}
|
||||
|
||||
// inverse fft (in-place)
|
||||
void ifft(CArray& x)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -0,0 +1,14 @@
|
|||
(defun fft (x)
|
||||
(if (<= (length x) 1) x
|
||||
(let*
|
||||
(
|
||||
(even (fft (loop for i from 0 below (length x) by 2 collect (nth i x))))
|
||||
(odd (fft (loop for i from 1 below (length x) by 2 collect (nth i x))))
|
||||
(aux (loop for k from 0 below (/ (length x) 2) collect (* (exp (/ (* (complex 0 -2) pi k ) (length x))) (nth k odd))))
|
||||
)
|
||||
(append (mapcar #'+ even aux) (mapcar #'- even aux))
|
||||
)
|
||||
)
|
||||
)
|
||||
|
||||
(mapcar (lambda (x) (format t "~a~&" x)) (fft '(1 1 1 1 0 0 0 0)))
|
||||
|
|
@ -29,9 +29,7 @@ contains
|
|||
|
||||
! combine
|
||||
do i=1,N/2
|
||||
|
||||
|
||||
t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),KIND=DP))*even(i)
|
||||
t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
|
||||
x(i) = odd(i) + t
|
||||
x(i+N/2) = odd(i) - t
|
||||
end do
|
||||
|
|
|
|||
|
|
@ -0,0 +1,22 @@
|
|||
package fft
|
||||
|
||||
import java.lang.Math.*
|
||||
|
||||
class Complex(val re: Double, val im: Double) {
|
||||
operator infix fun plus(x: Complex) = Complex(re + x.re, im + x.im)
|
||||
operator infix fun minus(x: Complex) = Complex(re - x.re, im - x.im)
|
||||
operator infix fun times(x: Double) = Complex(re * x, im * x)
|
||||
operator infix fun times(x: Complex) = Complex(re * x.re - im * x.im, re * x.im + im * x.re)
|
||||
operator infix fun div(x: Double) = Complex(re / x, im / x)
|
||||
val exp: Complex by lazy { Complex(cos(im), sin(im)) * (cosh(re) + sinh(re)) }
|
||||
|
||||
override fun toString() = when {
|
||||
b == "0.000" -> a
|
||||
a == "0.000" -> b + 'i'
|
||||
im > 0 -> a + " + " + b + 'i'
|
||||
else -> a + " - " + b + 'i'
|
||||
}
|
||||
|
||||
private final val a = "%1.3f".format(re)
|
||||
private final val b = "%1.3f".format(abs(im))
|
||||
}
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
package fft
|
||||
|
||||
object FFT {
|
||||
fun fft(a: Array<Complex>) = _fft(a, Complex(0.0, 2.0), 1.0)
|
||||
fun rfft(a: Array<Complex>) = _fft(a, Complex(0.0, -2.0), 2.0)
|
||||
|
||||
private fun _fft(a: Array<Complex>, direction: Complex, scalar: Double): Array<Complex> =
|
||||
if (a.size == 1)
|
||||
a
|
||||
else {
|
||||
val n = a.size
|
||||
require(n % 2 == 0, { "The Cooley-Tukey FFT algorithm only works when the length of the input is even." })
|
||||
|
||||
var (evens, odds) = Pair(emptyArray<Complex>(), emptyArray<Complex>())
|
||||
for (i in a.indices)
|
||||
if (i % 2 == 0) evens += a[i]
|
||||
else odds += a[i]
|
||||
evens = _fft(evens, direction, scalar)
|
||||
odds = _fft(odds, direction, scalar)
|
||||
|
||||
val pairs = (0 until n / 2).map {
|
||||
val offset = (direction * (java.lang.Math.PI * it / n)).exp * odds[it] / scalar
|
||||
val base = evens[it] / scalar
|
||||
Pair(base + offset, base - offset)
|
||||
}
|
||||
var (left, right) = Pair(emptyArray<Complex>(), emptyArray<Complex>())
|
||||
for ((l, r) in pairs) { left += l; right += r }
|
||||
left + right
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
package fft
|
||||
|
||||
fun Array<*>.println() = println(joinToString(prefix = "[", postfix = "]"))
|
||||
|
||||
fun main(args: Array<String>) {
|
||||
val data = arrayOf(Complex(1.0, 0.0), Complex(1.0, 0.0), Complex(1.0, 0.0), Complex(1.0, 0.0),
|
||||
Complex(0.0, 0.0), Complex(0.0, 2.0), Complex(0.0, 0.0), Complex(0.0, 0.0))
|
||||
|
||||
val a = FFT.fft(data)
|
||||
a.println()
|
||||
FFT.rfft(a).println()
|
||||
}
|
||||
|
|
@ -2,8 +2,8 @@ sub fft {
|
|||
return @_ if @_ == 1;
|
||||
my @evn = fft( @_[0, 2 ... *] );
|
||||
my @odd = fft( @_[1, 3 ... *] ) Z*
|
||||
(1, * * cis( 2 * pi / @_ ) ... *);
|
||||
return @evn »+« @odd, @evn »-« @odd;
|
||||
map &cis, (0, 2 * pi / @_ ... *);
|
||||
return flat @evn »+« @odd, @evn »-« @odd;
|
||||
}
|
||||
|
||||
my @seq = ^16;
|
||||
|
|
|
|||
|
|
@ -3,11 +3,11 @@ from cmath import exp, pi
|
|||
def fft(x):
|
||||
N = len(x)
|
||||
if N <= 1: return x
|
||||
even = fft2(x[0::2])
|
||||
odd = fft2(x[1::2])
|
||||
T= [exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]
|
||||
return [even[k] + T[k] for k in xrange(N/2)] + \
|
||||
[even[k] - T[k] for k in xrange(N/2)]
|
||||
even = fft(x[0::2])
|
||||
odd = fft(x[1::2])
|
||||
T= [exp(-2j*pi*k/N)*odd[k] for k in range(N//2)]
|
||||
return [even[k] + T[k] for k in range(N//2)] + \
|
||||
[even[k] - T[k] for k in range(N//2)]
|
||||
|
||||
print( ' '.join("%5.3f" % abs(f)
|
||||
for f in fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])) )
|
||||
|
|
|
|||
|
|
@ -1,81 +1,70 @@
|
|||
/*REXX pgm does a fast Fourier transform (FFT) on a set of complex nums.*/
|
||||
numeric digits length( pi() ) - 1 /*limited by PI function result. */
|
||||
arg data /*the ARG verb uppercases DATA */
|
||||
if data='' then data=1 1 1 1 0 /*No data? Then use the default.*/
|
||||
data=translate(data, 'J', "I") /*allow use of I as well as J */
|
||||
size=words(data); pad=left('',6) /*PAD: for indenting/padding SAYs*/
|
||||
do p=0 until 2**p>=size ; end /* # args exactly a power of 2? */
|
||||
do j=size+1 to 2**p;data=data 0; end /*add zeroes until a power of 2. */
|
||||
size=words(data); ph=p%2; call hdr /*┌─────────────────────────────┐*/
|
||||
/*│ Numbers in data can be in │*/
|
||||
do j=0 for size /*│ 7 formats: real │*/
|
||||
_=word(data,j+1) /*│ real,imag │*/
|
||||
parse var _ #.1.j ',' #.2.j /*│ ,imag │*/
|
||||
if right(#.1.j,1)=='J' then parse , /*│ nnnJ │*/
|
||||
var #.1.j #2.j "J" @.1.j /*│ nnnj │*/
|
||||
do m=1 for 2 /*omitted?*/ /*│ nnnI │*/
|
||||
#.m.j=word(#.m.j 0, 1) /*│ nnni │*/
|
||||
end /*m*/ /*└─────────────────────────────┘*/
|
||||
/*REXX pgm performs a fast Fourier transform (FFT) on a set of complex numbers*/
|
||||
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*/ /*║ Numbers in data can be in ║*/
|
||||
do j=0 for size /*║ seven formats: real ║*/
|
||||
_=translate(word(data,j+1), 'J', "I") /*║ real,imag ║*/
|
||||
parse var _ #.1.j '' $ 1 ',' #.2.j /*║ ,imag ║*/
|
||||
if $=='J' then parse var #.1.j #2.j , /*║ nnnJ ║*/
|
||||
"J" #.1.j /*║ nnnj ║*/
|
||||
do m=1 for 2; #.m.j=word(#.m.j 0,1) /*║ nnnI ║*/
|
||||
end /*m*/ /* [↑] ommited part?*/ /*║ nnni ║*/
|
||||
/*╚═════════════════════════════╝*/
|
||||
say pad " FFT in " center(j+1,7) pad fmt(#.1.j) fmt(#.2.j,'i')
|
||||
end /*j*/
|
||||
say
|
||||
tran=pi()*2/2**p; !.=0; hp=2**p%2; A=2**(p-ph); ptr=A; dbl=1
|
||||
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*/
|
||||
|
||||
say pad " FFT in " center(j+1,7) pad nice(#.1.j) nice(#.2.j,'i')
|
||||
end /*j*/
|
||||
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
|
||||
|
||||
say; say; tran=2*pi()/2**p; !.=0
|
||||
hp=2**p%2; counterA=2**(p-ph); pointer=counterA; doubler=1
|
||||
|
||||
do p-ph; halfpointer=pointer%2
|
||||
|
||||
do i=halfpointer by pointer to counterA-halfpointer
|
||||
_=i-halfpointer; !.i=!._+doubler
|
||||
end /*i*/
|
||||
|
||||
doubler=doubler*2; pointer=halfpointer
|
||||
end /*p-ph*/
|
||||
|
||||
do j=0 to 2**p%4; cmp.j=cos(j*tran); _m=hp-j; cmp._m=-cmp.j
|
||||
_p=hp+j; cmp._p=-cmp.j
|
||||
end /*j*/
|
||||
|
||||
counterB=2**ph
|
||||
|
||||
do i=0 for counterA; q=i *counterB
|
||||
do j=0 for counterB; h=q+j; _=!.j*counterB+!.i; if _<=h then iterate
|
||||
parse value #.1._ #.1.h #.2._ #.2.h with #.1.h #.1._ #.2.h #.2._
|
||||
end /*j*/ /* [↑] switch two sets of values*/
|
||||
do i=0 for A; q=i * B
|
||||
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._
|
||||
end /*j*/ /* [↑] swap two sets of values.*/
|
||||
end /*i*/
|
||||
|
||||
double=1; do p ; w=hp%double
|
||||
do k=0 for double; lb=w*k ; lh=lb+2**p%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 /*j*/
|
||||
end /*k*/
|
||||
double=double+double
|
||||
end /*p*/
|
||||
dbl=1; do p ; w=hp % dbl
|
||||
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
|
||||
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 /*j*/
|
||||
end /*k*/
|
||||
dbl=dbl+dbl
|
||||
end /*p*/
|
||||
call hdr
|
||||
do i=0 for size
|
||||
say pad " FFT out " center(i+1,7) pad nice(#.1.i) nice(#.2.i,'j')
|
||||
end /*i*/
|
||||
exit /*stick a fork in it, we're done.*/
|
||||
/*──────────────────────────────────HDR subroutine──────────────────────*/
|
||||
hdr: _='───data─── num real-part imaginary-part'; say pad _
|
||||
say pad translate(_, " "copies('═',256), " "xrange()); return
|
||||
/*──────────────────────────────────PI subroutine───────────────────────────────────*/
|
||||
pi: return , /*add more digs if NUMERIC DIGITS > 85. */
|
||||
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628
|
||||
/*──────────────────────────────────R2R subroutine──────────────────────*/
|
||||
r2r: return arg(1) // (2*pi()) /*reduce radians to unit circle. */
|
||||
/*──────────────────────────────────COS subroutine──────────────────────*/
|
||||
cos: procedure; parse 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
|
||||
/*──────────────────────────────────NICE subroutine─────────────────────*/
|
||||
nice: procedure; parse arg x,j /*makes complex nums look nicer. */
|
||||
numeric digits digits()%10; nz='1e-'digits() /*show ≈10% of 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)
|
||||
do i=0 for size
|
||||
say pad " FFT out " center(i+1,7) pad fmt(#.1.i) fmt(#.2.i,'j')
|
||||
end /*i*/ /*numbers are shown with 10 digs [↑] */
|
||||
exit /*stick a fork in it, we're all done. */
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
cos: procedure; parse arg x; q=r2r(x)**2; z=1; _=1; p=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 /*transforms complex numbers for looks.*/
|
||||
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,,.); if y>=0 then y=' 'y; return left(y||j, 12)
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
hdr: _='───data─── num real─part imaginary─part'; say pad _
|
||||
say pad translate(_, " "copies('═',256), " "xrange()); return
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
pi: return 3.141592653589793238462643383279502884197169399375105820974944592308
|
||||
/*────────────────────────────────────────────────────────────────────────────*/
|
||||
r2r: return arg(1) // (pi()*2) /*reduce the radians to a unit circle. */
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue