implement FFT;

include "short_fft.m";

include "draw.m";
draw : Draw;

include "math.m";
math : Math;

include "sys.m";
sys: Sys;

include "bufio.m";
bufio : Bufio;
Iobuf : import bufio;


init(nil: ref Draw->Context, nil: list of string) {
	sys = load Sys Sys->PATH;
	math = load Math Math->PATH;
	bufio = load Bufio Bufio->PATH;
}

Complexes.bit_reverse(c : self ref Complexes)
{
	tr, tj : real;
	i2 : int;
	i2 = (len c.r) >> 1;
 	j := 0;
	i : int;
	for (i = 0; i < (len c.r - 1); i++) {
		if (i < j) {
			tr = c.r[i];
			tj = c.j[i];
			c.r[i] = c.r[j];
			c.j[i] = c.j[j];
			c.r[j] = tr;
			c.j[j] = tj;
		}
		k := i2;
		while (k <= j) {
			if(k == 0) {
				sys->print("ACK WRONG NUMBER OF SAMPLES\nneeded %d got %d\n", 1 + i << 1, len c.r);return;
			}
			j -= k;
			k >>= 1;
		}
		j += k;
	}
}

Complexes.print(c : self ref Complexes, labelr, labeli : string)
{
	sys->print("wtf\n");
	for(k := 0; k < len c.r; k++)
		sys->print("k %d %s %f %s %f\n", k, labelr, c.r[k], labeli, c.j[k]);
}

Complexes.save(c : self ref Complexes, fd : ref Sys->FD, directory : string, frame : big, x, y : int) : ref Sys->FD
{
	if(fd == nil) {
		filename := sys->sprint("%s/%bd.%d-%d.fft", directory, frame, x, y);
		fd = sys->create(filename, sys->OWRITE, 0666);
		if(fd == nil) return nil;
	}

	b := array[len c.r * 8] of byte;
	math->export_real(b, c.r);
	sys->write(fd, b, len b);
	math->export_real(b, c.j);
	sys->write(fd, b, len b);
	
	return fd;
}

C_from_file(filename : string, num_channels, channel_no : int) : ref Complexes
{
	fd := sys->open(filename, sys->OREAD);
	if(fd == nil) return nil;
	(success, dir) := sys->fstat(fd);
	if(success < 0) return nil;
	size := int (dir.length >> 3); # 8 bytes per real

	sample_pairs := size / num_channels;
	samples := sample_pairs >> 1;

	sys->seek(fd, big (sample_pairs * channel_no), sys->SEEKSTART);

	r := array[samples] of real;
	j := array[samples] of real;

	b := array[8 * samples] of byte;
	sys->read(fd, b, len b);
	math->import_real(b, r);
	sys->read(fd, b, len b);
	math->import_real(b, j);

	c := ref Complexes(r, j);
	return c;
}

Complexes.fft(c : self ref Complexes, dir, m : int)
{
   n,i,i1,j,l,l1,l2 : int;
   c1,c2,t1,t2,u1,u2,z : real;

   # Calculate the number of points
   n =  1;
   for (i=0;i<m;i++) 
      n *= 2;

  # Do the bit reversal 

  c.bit_reverse();
   #  Compute the FFT 
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * c.r[i1] - u2 * c.j[i1];
            t2 = u1 * c.j[i1] + u2 * c.r[i1];
            c.r[i1] = c.r[i] - t1; 
            c.j[i1] = c.j[i] - t2;
            c.r[i] += t1;
            c.j[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = math->sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = math->sqrt((1.0 + c1) / 2.0);
   }

   #  Scaling for forward transform 
   if (dir == 1) {
      for (i=0;i<n;i++) {
         c.r[i] /= real n;
         c.j[i] /= real n;
      }
   }
   
}

C_from_ints(ints : array of int, maxvalue : real) : ref Complexes
{
	c := ref Complexes(array[len ints] of real, array[len ints] of real);
	for(i := 0; i < len ints; i++) {
		c.r[i] = real ints[i] / maxvalue;
		c.j[i] = 0.0;
	}
	return c;
}
C_from_reals(reals : array of real) : ref Complexes
{
	c := ref Complexes(reals, array[len reals] of real);
	for(i := 0; i < len reals; i++) {
		c.j[i] = 0.0;
	}
	return c;
}

Complexes.interpolate(c : self ref Complexes, k : ref Complexes, sc : real)
{
	if(len c.r != len k.r) {
		sys->print("MISMATCHED SAMPLES SIZES\n");
		return;
	}
	sk := 1.0 - sc;

	for(i := 0; i < len c.r; i++) {
		c.r[i] = sc * c.r[i] + sk * k.r[i];
		c.j[i] = sc * c.j[i] + sk * k.j[i];
	}
}