#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

float DoTest1(void);
void DoTest2(void);
void DoTest3(void);
void DoTest4(void);
void DoTest5(void);
void DoTest6(void);
void DoTest7(void);
void DoTest8(void);

int *ia;
float *a, *b;

int main(void)
{
clock_t tim;

	ia = (int *)malloc(32768 * sizeof(int));
	a = (float *)malloc(32768 * sizeof(float));
	b = (float *)malloc(32768 * sizeof(float));

	printf("Float matrix multiplication: ");
	tim = clock();
	DoTest1();
	printf("%d\n", clock() - tim);
	
	printf("Integer conditional move: ");
	tim = clock();
	DoTest2();
	printf("%d\n", clock() - tim);

	printf("Float conditional move: ");
	tim = clock();
	DoTest3();
	printf("%d\n", clock() - tim);

	printf("Lookup table accessing: ");
	tim = clock();
	DoTest4();
	printf("%d\n", clock() - tim);

	printf("Float biquad filter: ");
	tim = clock();
	DoTest5();
	printf("%d\n", clock() - tim);

	printf("Computing biquad coefficents: ");
	tim = clock();
	DoTest6();
	printf("%d\n", clock() - tim);

	printf("Int biquad filter: ");
	tim = clock();
	DoTest7();
	printf("%d\n", clock() - tim);

	printf("FFT: ");
	tim = clock();
	DoTest8();
	printf("%d\n", clock() - tim);

	free(b); free(a);
	free(ia);
	getchar();
	return(0);
}


float DoTest1(void)
{
float a[64][64], b[64][64], c[64][64];
float sum;
int i, j, k, count;

	for (count = 0; count < 1000; count++)
	{
		for (i = 0; i < 64; i++)
			for (j = 0; j < 64; j++)
				a[i][j] = b[i][j] = 1.0F;

		for (i = 0; i < 64; i++)
			for (j = 0; j < 64; j++)
			{
				sum = 0.0F;
				for (k = 0; k < 64; k++) sum += a[i][k] * b[k][j];
				c[i][j] = sum;
			}
	}
	return(c[5][6]);
}

void DoTest2(void)
{
int i, count;

	for (i = 0; i < 16384; i++)
		ia[i] = rand();
	
	for (count = 0; count < 10000; count++)
	{
		for (i = 0; i < 16384; i++)
			if (ia[i] >= 1)
				ia[i] = 0;
			else
				ia[i] = 1;
	}
}

void DoTest3(void)
{
int i, count;

	for (i = 0; i < 16384; i++)
		a[i] = (float)rand();
	
	for (count = 0; count < 10000; count++)
	{
		for (i = 0; i < 16384; i++)
			if (a[i] >= 1.0F)
				a[i] = 0.0F;
			else
				a[i] = 1.0F;
	}
}

void DoTest4(void)
{
int i, count;

	for (i = 0; i < 16384; i++)
		a[i] = ((float)(rand() % 16384)) / 16384.0F;

	for (count = 0; count < 1000; count++)
	{
		for (i = 0; i < 16384; i++)
			b[i] = a[(int)(a[i] * 16384.0F)];
	}
}

void DoTest5(void)
{
float a0, a1, a2, b1, b2;
float *ptIn, *ptOut;
register int i, count;

	for (i = 0; i < 16384; i++)
		a[i] = ((float)rand()) / 16384.0F;

	a2 = a0 = 0.0137F;
	a1 = 2.0F * a0;
	b1 = 2.0F * (a0 - 0.877F);
	b2 = 0.781F;
	for (count = 0; count < 5000; count++)
	{
		ptIn = a; ptOut = b;
		ptOut[0] = ptOut[1] = 0.0F;
		for (i = 2; i < 16384; i++)
			ptOut[i] = ptIn[i-0] * a0 + ptIn[i-1] * a1 + ptIn[i-2] * a2 -
								 ptOut[i-1] * b1 - ptOut[i-2] * b2;
	}
}


float Ga0, Ga1, Ga2, Gb1, Gb2;

void SetSPoly(float aNum, float bNum, float cNum, float aDen, float bDen, float cDen);
void SetSPoly(float aNum, float bNum, float cNum, float aDen, float bDen, float cDen)
{
float w = 2.0F * 44100.0F, w2;
float Norm;
float Num1, Num2, Den1, Den2;

	w2 = w * w;
	Num1 = w2 * aNum + cNum;
	Num2 = w * bNum;
	Den1 = w2 * aDen + cDen;
	Den2 = w * bDen;
	Norm = 1.0F / (Den1 + Den2);
	Ga0 = (Num1 + Num2) * Norm;
	Ga1 = 2.0F * (cNum - w2 * aNum) * Norm;
	Ga2 = (Num1 - Num2) * Norm;
	Gb1 = 2.0F * (cDen - w2 * aDen) * Norm;
	Gb2 = (Den1 - Den2) * Norm;
}

void DoTest6(void)
{
float f;
int i, count;

	for (count = 0; count < 50000; count++)
		for (i = 0, f = 100.0F; i < 500; i++, f *= 1.01F)
			SetSPoly(1.0F, 0.0F, 0.0F, 1.0F, f, f * f);
}


void DoTest7(void)
{
int a0, a1, a2, b1, b2;
int *ptIn, *ptOut;
register int i, count;

	for (i = 0; i < 16384; i++)
		a[i] = ((float)rand()) / 16384.0F;

	a2 = a0 = (int)(0.0137F * (1 << 12));
	a1 = 2 * a0;
	b1 = 2 * (a0 - (int)(0.877F * (1 << 12)));
	b2 = (int)(0.781F * (1 << 12));
	for (count = 0; count < 5000; count++)
	{
		ptIn = (int *)a; ptOut = (int *)b;
		ptOut[0] = ptOut[1] = 0;
		for (i = 2; i < 16384; i++)
			ptOut[i] = (ptIn[i-0] * a0 + ptIn[i-1] * a1 + ptIn[i-2] * a2 -
								 ptOut[i-1] * b1 - ptOut[i-2] * b2) >> 12;
	}
}



/*
	The code of fft is by Don Cross <dcross@intersrv.com>
	http://www.intersrv.com/~dcross/fft.html
*/

unsigned NumberOfBitsNeeded(unsigned PowerOfTwo)
{
	for (unsigned i = 0;; i++)
		if (PowerOfTwo & (1 << i))
			 return i;
}

unsigned ReverseBits (unsigned index, unsigned NumBits)
{
unsigned i, rev;

	for (i = rev = 0; i < NumBits; i++)
  {
    rev = (rev << 1) | (index & 1);
    index >>= 1;
	}
  return rev;
}

void fft_float ( unsigned NumSamples,
                 int    InverseTransform,
                 float  *RealIn,
                 float  *ImagIn,
                 float  *RealOut,
                 float  *ImagOut )
{
unsigned NumBits;    /* Number of bits needed to store indices */
unsigned i, j, k, n;
unsigned BlockSize, BlockEnd;

float angle_numerator = 2.0F * 3.141592F;
float delta_angle;
float alpha, beta;  /* used in recurrence relation */
float delta_ar;
float tr, ti;     /* temp real, temp imaginary */
float ar, ai;     /* angle vector real, angle vector imaginary */

	if (InverseTransform)
		angle_numerator = -angle_numerator;
	NumBits = NumberOfBitsNeeded ( NumSamples );

	for ( i=0; i < NumSamples; i++ )
	{
		j = ReverseBits ( i, NumBits );
		RealOut[j] = RealIn[i];
		ImagOut[j] = ImagIn ? ImagIn[i] : 0.0F;
	}

	BlockEnd = 1;
	for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
	{
		delta_angle = angle_numerator / (float)BlockSize;
		alpha = (float)sin(0.5 * delta_angle);
		alpha = 2.0F * alpha * alpha;
		beta = (float)sin (delta_angle);

		for ( i=0; i < NumSamples; i += BlockSize )
		{
			ar = 1.0;   /* cos(0) */
			ai = 0.0;   /* sin(0) */
			 for ( j=i, n=0; n < BlockEnd; j++, n++ )
			 {
					k = j + BlockEnd;
					tr = ar*RealOut[k] - ai*ImagOut[k];
					ti = ar*ImagOut[k] + ai*RealOut[k];

					RealOut[k] = RealOut[j] - tr;
					ImagOut[k] = ImagOut[j] - ti;

					RealOut[j] += tr;
					ImagOut[j] += ti;

					delta_ar = alpha*ar + beta*ai;
					ai -= (alpha*ai - beta*ar);
					ar -= delta_ar;
			 }
		}

		BlockEnd = BlockSize;
	}

	if ( InverseTransform )
	{
		float denom = (float)NumSamples;
		for ( i=0; i < NumSamples; i++ )
		{
			 RealOut[i] /= denom;
			 ImagOut[i] /= denom;
		}
	}
}


void DoTest8(void)
{
register int i, count;

	for (i = 0; i < 16384; i++)
		a[i] = ((float)rand()) / 16384.0F;

	for (count = 0; count < 500; count++)
		fft_float(8192, 0, a, a + 8192, b, b + 8192);
}