#define	SAMPLES		256
#define	POWER		8	/* 2 to the 8th power is 256 */


/*
 *	Fixed constants should be used in the macros for speed.
 *
 *	A cosine wave leads a sine wave by 90 degrees, so offset into
 *	the lookup table 1/4 the way into it (256 / 4 = 64).  Using
 *	modulo 256, lookups will wrap around to zero for numbers greater
 *	than 255.  (cosine(200) = Sine[264 % 256] = Sine[8]).
 */


#define	sine(x)		Sine[(x)]
#define cosine(x)	Sine[((x) + 64) % 256]

/* Globals, Forward declarations */

double	Real[SAMPLES], Imag[SAMPLES], Maxn, magnitude();
int	GraphDriver = DETECT, GraphMode, Primary, Cursor, Length;
int	Bottom, Left, PrintChar(), PrintScreen(), getkey();
void	*Save, quit(), beep(), build_window(), commands(), bee_bop();
void	scale(), fft(), max_amp(), display();

/*
 *	Bit Reverse Table for size 256
 *	Lookup saves 20 seconds in Turbo-C over the pow() function.
 *
 *	Br_table[x] = x inverted (eg. 00000001 flipped to 10000000)
 */

unsigned char permute[] = {
	0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0, 0x10, 0x90, 0x50, 0xd0,
	0x30, 0xb0, 0x70, 0xf0, 0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8,
	0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8, 0x04, 0x84, 0x44, 0xc4,
	0x24, 0xa4, 0x64, 0xe4, 0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4,
	0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec, 0x1c, 0x9c, 0x5c, 0xdc,
	0x3c, 0xbc, 0x7c, 0xfc, 0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2,
	0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2, 0x0a, 0x8a, 0x4a, 0xca,
	0x2a, 0xaa, 0x6a, 0xea, 0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa,
	0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6, 0x16, 0x96, 0x56, 0xd6,
	0x36, 0xb6, 0x76, 0xf6, 0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee,
	0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe, 0x01, 0x81, 0x41, 0xc1,
	0x21, 0xa1, 0x61, 0xe1, 0x11, 0x91, 0x51, 0xd1, 0x31, 0xb1, 0x71, 0xf1,
	0x09, 0x89, 0x49, 0xc9, 0x29, 0xa9, 0x69, 0xe9, 0x19, 0x99, 0x59, 0xd9,
	0x39, 0xb9, 0x79, 0xf9, 0x05, 0x85, 0x45, 0xc5, 0x25, 0xa5, 0x65, 0xe5,
	0x15, 0x95, 0x55, 0xd5, 0x35, 0xb5, 0x75, 0xf5, 0x0d, 0x8d, 0x4d, 0xcd,
	0x2d, 0xad, 0x6d, 0xed, 0x1d, 0x9d, 0x5d, 0xdd, 0x3d, 0xbd, 0x7d, 0xfd,
	0x03, 0x83, 0x43, 0xc3, 0x23, 0xa3, 0x63, 0xe3, 0x13, 0x93, 0x53, 0xd3,
	0x33, 0xb3, 0x73, 0xf3, 0x0b, 0x8b, 0x4b, 0xcb, 0x2b, 0xab, 0x6b, 0xeb,
	0x1b, 0x9b, 0x5b, 0xdb, 0x3b, 0xbb, 0x7b, 0xfb, 0x07, 0x87, 0x47, 0xc7,
	0x27, 0xa7, 0x67, 0xe7, 0x17, 0x97, 0x57, 0xd7, 0x37, 0xb7, 0x77, 0xf7,
	0x0f, 0x8f, 0x4f, 0xcf, 0x2f, 0xaf, 0x6f, 0xef, 0x1f, 0x9f, 0x5f, 0xdf,
	0x3f, 0xbf, 0x7f, 0xff
};

/*
 *	Sine/Cosine Table for size 256, Lookup saves 9 seconds in Turbo-C
 *
 *	Sine[n] = sin(x), x = x + (2Pi / 256)
 */

float	Sine[] = {
	0.000000, 0.024541, 0.049068, 0.073565,	0.098017, 0.122411,
	0.146730, 0.170962, 0.195090, 0.219101, 0.242980, 0.266713,
	0.290285, 0.313682, 0.336890, 0.359895, 0.382683, 0.405241,
	0.427555, 0.449611, 0.471397, 0.492898, 0.514103, 0.534998,
	0.555570, 0.575808, 0.595699, 0.615232,	0.634393, 0.653173,
	0.671559, 0.689541, 0.707107, 0.724247, 0.740951, 0.757209,
	0.773010, 0.788346, 0.803208, 0.817585, 0.831470, 0.844854,
	0.857729, 0.870087, 0.881921, 0.893224, 0.903989, 0.914210,
	0.923880, 0.932993, 0.941544, 0.949528,	0.956940, 0.963776,
	0.970031, 0.975702, 0.980785, 0.985278, 0.989177, 0.992480,
	0.995185, 0.997290, 0.998795, 0.999699,	1.000000, 0.999699,
	0.998795, 0.997290, 0.995185, 0.992480, 0.989177, 0.985278,
	0.980785, 0.975702, 0.970031, 0.963776,	0.956940, 0.949528,
	0.941544, 0.932993, 0.923880, 0.914210, 0.903989, 0.893224,
	0.881921, 0.870087, 0.857729, 0.844854,	0.831470, 0.817585,
	0.803208, 0.788346, 0.773010, 0.757209, 0.740951, 0.724247,
	0.707107, 0.689541, 0.671559, 0.653173,	0.634393, 0.615232,
	0.595699, 0.575808, 0.555570, 0.534998, 0.514103, 0.492898,
	0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895,
	0.336890, 0.313682, 0.290285, 0.266713, 0.242980, 0.219101,
	0.195090, 0.170962, 0.146730, 0.122411,	0.098017, 0.073565,
	0.049068, 0.024541, 0.000000, -0.024541, -0.049068, -0.073565,
	-0.098017, -0.122411, -0.146730, -0.170962, -0.195090, -0.219101,
	-0.242980, -0.266713, -0.290285, -0.313682, -0.336890, -0.359895,
	-0.382683, -0.405241, -0.427555, -0.449611, -0.471397, -0.492898,
	-0.514103, -0.534998, -0.555570, -0.575808, -0.595699, -0.615232,
	-0.634393, -0.653173, -0.671559, -0.689541, -0.707107, -0.724247,
	-0.740951, -0.757209, -0.773010, -0.788346, -0.803208, -0.817585,
	-0.831470, -0.844854, -0.857729, -0.870087, -0.881921, -0.893224,
	-0.903989, -0.914210, -0.923880, -0.932993, -0.941544, -0.949528,
	-0.956940, -0.963776, -0.970031, -0.975702, -0.980785, -0.985278,
	-0.989177, -0.992480, -0.995185, -0.997290, -0.998795, -0.999699,
	-1.000000, -0.999699, -0.998795, -0.997290, -0.995185, -0.992480,
	-0.989177, -0.985278, -0.980785, -0.975702, -0.970031, -0.963776,
	-0.956940, -0.949528, -0.941544, -0.932993, -0.923880, -0.914210,
	-0.903989, -0.893224, -0.881921, -0.870087, -0.857729, -0.844854,
	-0.831470, -0.817585, -0.803208, -0.788346, -0.773010, -0.757209,
	-0.740951, -0.724247, -0.707107, -0.689541, -0.671559, -0.653173,
	-0.634393, -0.615232, -0.595699, -0.575808, -0.555570, -0.534998,
	-0.514103, -0.492898, -0.471397, -0.449611, -0.427555, -0.405241,
	-0.382683, -0.359895, -0.336890, -0.313682, -0.290285, -0.266713,
	-0.242980, -0.219101, -0.195090, -0.170962, -0.146730, -0.122411,
	-0.098017, -0.073565, -0.049068, -0.024541
};



void scale()
{
	register int	loop;

	for (loop = 0; loop < SAMPLES; loop++)  {
		Real[loop] /= SAMPLES;
		Imag[loop] /= SAMPLES;
	}
}


void fft()
{
	register int	loop, loop1, loop2;
	unsigned	i1;			/* going to right shift this */
	int		i2, i3, i4, y;
	double		a1, a2, b1, b2, z1, z2;

	i1 = SAMPLES >> 1;
	i2 = 1;

	/* perform the butterfly's */

	for (loop = 0; loop < POWER; loop++)  {
		i3 = 0;
		i4 = i1;

		for (loop1 = 0; loop1 < i2; loop1++)  {
			y = permute(i3 / (int)i1);
			z1 =  cosine(y);
			z2 = -sine(y);

			for (loop2 = i3; loop2 < i4; loop2++)  {

				a1 = Real[loop2];
				a2 = Imag[loop2];

				b1  = z1*Real[loop2+i1] - z2*Imag[loop2+i1];
				b2  = z2*Real[loop2+i1] + z1*Imag[loop2+i1];

				Real[loop2] = a1 + b1;
				Imag[loop2] = a2 + b2;

				Real[loop2+i1] = a1 - b1;
				Imag[loop2+i1] = a2 - b2;
			}

			i3 += (i1 << 1);
			i4 += (i1 << 1);
		}

		i1 >>= 1;
		i2 <<= 1;
	}
}




/* Find maximum amplitude */
void max_amp()
{
	register int	loop;
	double		mag;

	Maxn = 0.0;
	for (loop = 0; loop < SAMPLES; loop++)  {
		if ((mag = magnitude(loop)) > Maxn)
			Maxn = mag;
	}
}

/*
 *	Calculate Power Magnitude
 */

double magnitude(n)
int	n;
{
	n = permute(n);
	return hypot(Real[n], Imag[n]);
}