// FFT.cpp : ƒCƒ“ƒvƒŠƒ�ƒ“ƒe�[ƒVƒ‡ƒ“ ƒtƒ@ƒCƒ‹
//

#import "FFT.h"

static void rdft(int n, int isgn, float *a, int *ip, float *w);
static void makewt(int nw,	int *ip, float *w);
static void makect(int nc, int *ip, float *c);
static void bitrv2(int n, int *ip, float *a);
static void cftfsub(int n,	float *a, const float *w);
static void cftbsub(int n,	float *a, const float *w);
static void cft1st(int n, float *a, const float *w);
static void cftmdl(int n, int l, float *a, const float *w);
static void rftfsub(int n, float *a, int nc, const float *c);
static void rftbsub(int n, float *a, int nc, const float *c);

CRFFT::CRFFT()
{
	last_n = 0;
	a = nullptr;
	ip = NULL;
	w = NULL;
}

CRFFT::~CRFFT()
{
	alloc_work(0);
}

void CRFFT::fft(int n, float *fp)
{
	fft_sub(n, fp, 1);
}

void CRFFT::ifft(int n, float *fp)
{
	fft_sub(n, fp, -1);
}

void CRFFT::fft_sub(int n, float *fp, int isgn)
{
	int i;

	alloc_work(n);

	if (n == 0)
		return;

	for (i = 0; i < n; i++)
		a[i] = fp[i];

	rdft(n, isgn, a, ip, w);

	for (i = 0; i < n; i++)
		fp[i] = a[i];
}

void CRFFT::alloc_work(int n)
{
	if (last_n != n) {
		if (a != NULL) {
			delete [] a;
			a = NULL;
		}

		if (ip != NULL) {
			delete [] ip;
			ip = NULL;
		}

		if (w != NULL) {
			delete [] w;
			w = NULL;
		}

		if (n != 0) {
			a = new float[n * 2];
			ip = new int[(int)sqrt((float)n) + 3];
			ip[0] = 0;
			w = new float[n / 2];
		}

		last_n = n;
	}
}

void rdft(int n, int isgn, float *a, int *ip, float *w)
{
    int nw, nc;
    float xi;

    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > (nc << 2)) {
        nc = n >> 2;
        makect(nc, ip, w + nw);
    }
    if (isgn >= 0) {
        if (n > 4) {
            bitrv2(n, ip + 2, a);
            cftfsub(n, a, w);
            rftfsub(n, a, nc, w + nw);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
        xi = a[0] - a[1];
        a[0] += a[1];
        a[1] = xi;
    } else {
        a[1] = 0.5 * (a[0] - a[1]);
        a[0] -= a[1];
        if (n > 4) {
            rftbsub(n, a, nc, w + nw);
            bitrv2(n, ip + 2, a);
            cftbsub(n, a, w);
        } else if (n == 4) {
            cftfsub(n, a, w);
        }
    }
}

/* --------	initializing routines -------- */


void makewt(int nw,	int *ip, float *w)
{
	int j, nwh;
	float delta, x, y;

	ip[0] =	nw;
	ip[1] =	1;
	if (nw > 2)	{
		nwh	= nw >>	1;
		delta =	atan(1.0) /	nwh;
		w[0] = 1;
		w[1] = 0;
		w[nwh] = cos(delta * nwh);
		w[nwh +	1] = w[nwh];
		if (nwh	> 2) {
			for	(j = 2;	j <	nwh; j += 2) {
				x =	cos(delta *	j);
				y =	sin(delta *	j);
				w[j] = x;
				w[j	+ 1] = y;
				w[nw - j] =	y;
				w[nw - j + 1] =	x;
			}
			bitrv2(nw, ip +	2, w);
		}
	}
}

void makect(int nc, int *ip, float *c)
{
    int j, nch;
    float delta;

    ip[1] = nc;
    if (nc > 1) {
        nch = nc >> 1;
        delta = atan(1.0) / nch;
        c[0] = cos(delta * nch);
        c[nch] = 0.5 * c[0];
        for (j = 1; j < nch; j++) {
            c[j] = 0.5 * cos(delta * j);
            c[nc - j] = 0.5 * sin(delta * j);
        }
    }
}


/* --------	child routines --------	*/


void bitrv2(int n, int *ip,	float *a)
{
	int j, j1, k, k1, l, m,	m2;
	float xr, xi, yr, yi;

	ip[0] =	0;
	l =	n;
	m =	1;
	while ((m << 3)	< l) {
		l >>= 1;
		for	(j = 0;	j <	m; j++)	{
			ip[m + j] =	ip[j] +	l;
		}
		m <<= 1;
	}
	m2 = 2 * m;
	if ((m << 3) ==	l) {
		for	(k = 0;	k <	m; k++)	{
			for	(j = 0;	j <	k; j++)	{
				j1 = 2 * j + ip[k];
				k1 = 2 * k + ip[j];
				xr = a[j1];
				xi = a[j1 +	1];
				yr = a[k1];
				yi = a[k1 +	1];
				a[j1] =	yr;
				a[j1 + 1] =	yi;
				a[k1] =	xr;
				a[k1 + 1] =	xi;
				j1 += m2;
				k1 += 2	* m2;
				xr = a[j1];
				xi = a[j1 +	1];
				yr = a[k1];
				yi = a[k1 +	1];
				a[j1] =	yr;
				a[j1 + 1] =	yi;
				a[k1] =	xr;
				a[k1 + 1] =	xi;
				j1 += m2;
				k1 -= m2;
				xr = a[j1];
				xi = a[j1 +	1];
				yr = a[k1];
				yi = a[k1 +	1];
				a[j1] =	yr;
				a[j1 + 1] =	yi;
				a[k1] =	xr;
				a[k1 + 1] =	xi;
				j1 += m2;
				k1 += 2	* m2;
				xr = a[j1];
				xi = a[j1 +	1];
				yr = a[k1];
				yi = a[k1 +	1];
				a[j1] =	yr;
				a[j1 + 1] =	yi;
				a[k1] =	xr;
				a[k1 + 1] =	xi;
			}
			j1 = 2 * k + m2	+ ip[k];
			k1 = j1	+ m2;
			xr = a[j1];
			xi = a[j1 +	1];
			yr = a[k1];
			yi = a[k1 +	1];
			a[j1] =	yr;
			a[j1 + 1] =	yi;
			a[k1] =	xr;
			a[k1 + 1] =	xi;
		}
	} else {
		for	(k = 1;	k <	m; k++)	{
			for	(j = 0;	j <	k; j++)	{
				j1 = 2 * j + ip[k];
				k1 = 2 * k + ip[j];
				xr = a[j1];
				xi = a[j1 +	1];
				yr = a[k1];
				yi = a[k1 +	1];
				a[j1] =	yr;
				a[j1 + 1] =	yi;
				a[k1] =	xr;
				a[k1 + 1] =	xi;
				j1 += m2;
				k1 += m2;
				xr = a[j1];
				xi = a[j1 +	1];
				yr = a[k1];
				yi = a[k1 +	1];
				a[j1] =	yr;
				a[j1 + 1] =	yi;
				a[k1] =	xr;
				a[k1 + 1] =	xi;
			}
		}
	}
}



void cftfsub(int n,	float *a, const float *w)
{
	int j, j1, j2, j3, l;
	float x0r,	x0i, x1r, x1i, x2r,	x2i, x3r, x3i;

	l =	2;
	if (n >	8) {
		cft1st(n, a, w);
		l =	8;
		while ((l << 2)	< n) {
			cftmdl(n, l, a,	w);
			l <<= 2;
		}
	}
	if ((l << 2) ==	n) {
		for	(j = 0;	j <	l; j +=	2) {
			j1 = j + l;
			j2 = j1	+ l;
			j3 = j2	+ l;
			x0r	= a[j] + a[j1];
			x0i	= a[j +	1] + a[j1 +	1];
			x1r	= a[j] - a[j1];
			x1i	= a[j +	1] - a[j1 +	1];
			x2r	= a[j2]	+ a[j3];
			x2i	= a[j2 + 1]	+ a[j3 + 1];
			x3r	= a[j2]	- a[j3];
			x3i	= a[j2 + 1]	- a[j3 + 1];
			a[j] = x0r + x2r;
			a[j	+ 1] = x0i + x2i;
			a[j2] =	x0r	- x2r;
			a[j2 + 1] =	x0i	- x2i;
			a[j1] =	x1r	- x3i;
			a[j1 + 1] =	x1i	+ x3r;
			a[j3] =	x1r	+ x3i;
			a[j3 + 1] =	x1i	- x3r;
		}
	} else {
		for	(j = 0;	j <	l; j +=	2) {
			j1 = j + l;
			x0r	= a[j] - a[j1];
			x0i	= a[j +	1] - a[j1 +	1];
			a[j] +=	a[j1];
			a[j	+ 1] +=	a[j1 + 1];
			a[j1] =	x0r;
			a[j1 + 1] =	x0i;
		}
	}
}


void cftbsub(int n,	float *a, const float *w)
{
	int j, j1, j2, j3, l;
	float x0r,	x0i, x1r, x1i, x2r,	x2i, x3r, x3i;

	l =	2;
	if (n >	8) {
		cft1st(n, a, w);
		l =	8;
		while ((l << 2)	< n) {
			cftmdl(n, l, a,	w);
			l <<= 2;
		}
	}
	if ((l << 2) ==	n) {
		for	(j = 0;	j <	l; j +=	2) {
			j1 = j + l;
			j2 = j1	+ l;
			j3 = j2	+ l;
			x0r	= a[j] + a[j1];
			x0i	= -a[j + 1]	- a[j1 + 1];
			x1r	= a[j] - a[j1];
			x1i	= -a[j + 1]	+ a[j1 + 1];
			x2r	= a[j2]	+ a[j3];
			x2i	= a[j2 + 1]	+ a[j3 + 1];
			x3r	= a[j2]	- a[j3];
			x3i	= a[j2 + 1]	- a[j3 + 1];
			a[j] = x0r + x2r;
			a[j	+ 1] = x0i - x2i;
			a[j2] =	x0r	- x2r;
			a[j2 + 1] =	x0i	+ x2i;
			a[j1] =	x1r	- x3i;
			a[j1 + 1] =	x1i	- x3r;
			a[j3] =	x1r	+ x3i;
			a[j3 + 1] =	x1i	+ x3r;
		}
	} else {
		for	(j = 0;	j <	l; j +=	2) {
			j1 = j + l;
			x0r	= a[j] - a[j1];
			x0i	= -a[j + 1]	+ a[j1 + 1];
			a[j] +=	a[j1];
			a[j	+ 1] = -a[j	+ 1] - a[j1	+ 1];
			a[j1] =	x0r;
			a[j1 + 1] =	x0i;
		}
	}
}


void cft1st(int n, float *a, const float *w)
{
	int j, k1, k2;
	float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
	float x0r,	x0i, x1r, x1i, x2r,	x2i, x3r, x3i;

	x0r	= a[0] + a[2];
	x0i	= a[1] + a[3];
	x1r	= a[0] - a[2];
	x1i	= a[1] - a[3];
	x2r	= a[4] + a[6];
	x2i	= a[5] + a[7];
	x3r	= a[4] - a[6];
	x3i	= a[5] - a[7];
	a[0] = x0r + x2r;
	a[1] = x0i + x2i;
	a[4] = x0r - x2r;
	a[5] = x0i - x2i;
	a[2] = x1r - x3i;
	a[3] = x1i + x3r;
	a[6] = x1r + x3i;
	a[7] = x1i - x3r;
	wk1r = w[2];
	x0r	= a[8] + a[10];
	x0i	= a[9] + a[11];
	x1r	= a[8] - a[10];
	x1i	= a[9] - a[11];
	x2r	= a[12]	+ a[14];
	x2i	= a[13]	+ a[15];
	x3r	= a[12]	- a[14];
	x3i	= a[13]	- a[15];
	a[8] = x0r + x2r;
	a[9] = x0i + x2i;
	a[12] =	x2i	- x0i;
	a[13] =	x0r	- x2r;
	x0r	= x1r -	x3i;
	x0i	= x1i +	x3r;
	a[10] =	wk1r * (x0r	- x0i);
	a[11] =	wk1r * (x0r	+ x0i);
	x0r	= x3i +	x1r;
	x0i	= x3r -	x1i;
	a[14] =	wk1r * (x0i	- x0r);
	a[15] =	wk1r * (x0i	+ x0r);
	k1 = 0;
	for	(j = 16; j < n;	j += 16) {
		k1 += 2;
		k2 = 2 * k1;
		wk2r = w[k1];
		wk2i = w[k1	+ 1];
		wk1r = w[k2];
		wk1i = w[k2	+ 1];
		wk3r = wk1r	- 2	* wk2i * wk1i;
		wk3i = 2 * wk2i	* wk1r - wk1i;
		x0r	= a[j] + a[j + 2];
		x0i	= a[j +	1] + a[j + 3];
		x1r	= a[j] - a[j + 2];
		x1i	= a[j +	1] - a[j + 3];
		x2r	= a[j +	4] + a[j + 6];
		x2i	= a[j +	5] + a[j + 7];
		x3r	= a[j +	4] - a[j + 6];
		x3i	= a[j +	5] - a[j + 7];
		a[j] = x0r + x2r;
		a[j	+ 1] = x0i + x2i;
		x0r	-= x2r;
		x0i	-= x2i;
		a[j	+ 4] = wk2r	* x0r -	wk2i * x0i;
		a[j	+ 5] = wk2r	* x0i +	wk2i * x0r;
		x0r	= x1r -	x3i;
		x0i	= x1i +	x3r;
		a[j	+ 2] = wk1r	* x0r -	wk1i * x0i;
		a[j	+ 3] = wk1r	* x0i +	wk1i * x0r;
		x0r	= x1r +	x3i;
		x0i	= x1i -	x3r;
		a[j	+ 6] = wk3r	* x0r -	wk3i * x0i;
		a[j	+ 7] = wk3r	* x0i +	wk3i * x0r;
		wk1r = w[k2	+ 2];
		wk1i = w[k2	+ 3];
		wk3r = wk1r	- 2	* wk2r * wk1i;
		wk3i = 2 * wk2r	* wk1r - wk1i;
		x0r	= a[j +	8] + a[j + 10];
		x0i	= a[j +	9] + a[j + 11];
		x1r	= a[j +	8] - a[j + 10];
		x1i	= a[j +	9] - a[j + 11];
		x2r	= a[j +	12]	+ a[j +	14];
		x2i	= a[j +	13]	+ a[j +	15];
		x3r	= a[j +	12]	- a[j +	14];
		x3i	= a[j +	13]	- a[j +	15];
		a[j	+ 8] = x0r + x2r;
		a[j	+ 9] = x0i + x2i;
		x0r	-= x2r;
		x0i	-= x2i;
		a[j	+ 12] =	-wk2i *	x0r	- wk2r * x0i;
		a[j	+ 13] =	-wk2i *	x0i	+ wk2r * x0r;
		x0r	= x1r -	x3i;
		x0i	= x1i +	x3r;
		a[j	+ 10] =	wk1r * x0r - wk1i *	x0i;
		a[j	+ 11] =	wk1r * x0i + wk1i *	x0r;
		x0r	= x1r +	x3i;
		x0i	= x1i -	x3r;
		a[j	+ 14] =	wk3r * x0r - wk3i *	x0i;
		a[j	+ 15] =	wk3r * x0i + wk3i *	x0r;
	}
}


void cftmdl(int n, int l, float *a, const float *w)
{
	int j, j1, j2, j3, k, k1, k2, m, m2;
	float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
	float x0r,	x0i, x1r, x1i, x2r,	x2i, x3r, x3i;

	m =	l << 2;
	for	(j = 0;	j <	l; j +=	2) {
		j1 = j + l;
		j2 = j1	+ l;
		j3 = j2	+ l;
		x0r	= a[j] + a[j1];
		x0i	= a[j +	1] + a[j1 +	1];
		x1r	= a[j] - a[j1];
		x1i	= a[j +	1] - a[j1 +	1];
		x2r	= a[j2]	+ a[j3];
		x2i	= a[j2 + 1]	+ a[j3 + 1];
		x3r	= a[j2]	- a[j3];
		x3i	= a[j2 + 1]	- a[j3 + 1];
		a[j] = x0r + x2r;
		a[j	+ 1] = x0i + x2i;
		a[j2] =	x0r	- x2r;
		a[j2 + 1] =	x0i	- x2i;
		a[j1] =	x1r	- x3i;
		a[j1 + 1] =	x1i	+ x3r;
		a[j3] =	x1r	+ x3i;
		a[j3 + 1] =	x1i	- x3r;
	}
	wk1r = w[2];
	for	(j = m;	j <	l +	m; j +=	2) {
		j1 = j + l;
		j2 = j1	+ l;
		j3 = j2	+ l;
		x0r	= a[j] + a[j1];
		x0i	= a[j +	1] + a[j1 +	1];
		x1r	= a[j] - a[j1];
		x1i	= a[j +	1] - a[j1 +	1];
		x2r	= a[j2]	+ a[j3];
		x2i	= a[j2 + 1]	+ a[j3 + 1];
		x3r	= a[j2]	- a[j3];
		x3i	= a[j2 + 1]	- a[j3 + 1];
		a[j] = x0r + x2r;
		a[j	+ 1] = x0i + x2i;
		a[j2] =	x2i	- x0i;
		a[j2 + 1] =	x0r	- x2r;
		x0r	= x1r -	x3i;
		x0i	= x1i +	x3r;
		a[j1] =	wk1r * (x0r	- x0i);
		a[j1 + 1] =	wk1r * (x0r	+ x0i);
		x0r	= x3i +	x1r;
		x0i	= x3r -	x1i;
		a[j3] =	wk1r * (x0i	- x0r);
		a[j3 + 1] =	wk1r * (x0i	+ x0r);
	}
	k1 = 0;
	m2 = 2 * m;
	for	(k = m2; k < n;	k += m2) {
		k1 += 2;
		k2 = 2 * k1;
		wk2r = w[k1];
		wk2i = w[k1	+ 1];
		wk1r = w[k2];
		wk1i = w[k2	+ 1];
		wk3r = wk1r	- 2	* wk2i * wk1i;
		wk3i = 2 * wk2i	* wk1r - wk1i;
		for	(j = k;	j <	l +	k; j +=	2) {
			j1 = j + l;
			j2 = j1	+ l;
			j3 = j2	+ l;
			x0r	= a[j] + a[j1];
			x0i	= a[j +	1] + a[j1 +	1];
			x1r	= a[j] - a[j1];
			x1i	= a[j +	1] - a[j1 +	1];
			x2r	= a[j2]	+ a[j3];
			x2i	= a[j2 + 1]	+ a[j3 + 1];
			x3r	= a[j2]	- a[j3];
			x3i	= a[j2 + 1]	- a[j3 + 1];
			a[j] = x0r + x2r;
			a[j	+ 1] = x0i + x2i;
			x0r	-= x2r;
			x0i	-= x2i;
			a[j2] =	wk2r * x0r - wk2i *	x0i;
			a[j2 + 1] =	wk2r * x0i + wk2i *	x0r;
			x0r	= x1r -	x3i;
			x0i	= x1i +	x3r;
			a[j1] =	wk1r * x0r - wk1i *	x0i;
			a[j1 + 1] =	wk1r * x0i + wk1i *	x0r;
			x0r	= x1r +	x3i;
			x0i	= x1i -	x3r;
			a[j3] =	wk3r * x0r - wk3i *	x0i;
			a[j3 + 1] =	wk3r * x0i + wk3i *	x0r;
		}
		wk1r = w[k2	+ 2];
		wk1i = w[k2	+ 3];
		wk3r = wk1r	- 2	* wk2r * wk1i;
		wk3i = 2 * wk2r	* wk1r - wk1i;
		for	(j = k + m;	j <	l +	(k + m); j += 2) {
			j1 = j + l;
			j2 = j1	+ l;
			j3 = j2	+ l;
			x0r	= a[j] + a[j1];
			x0i	= a[j +	1] + a[j1 +	1];
			x1r	= a[j] - a[j1];
			x1i	= a[j +	1] - a[j1 +	1];
			x2r	= a[j2]	+ a[j3];
			x2i	= a[j2 + 1]	+ a[j3 + 1];
			x3r	= a[j2]	- a[j3];
			x3i	= a[j2 + 1]	- a[j3 + 1];
			a[j] = x0r + x2r;
			a[j	+ 1] = x0i + x2i;
			x0r	-= x2r;
			x0i	-= x2i;
			a[j2] =	-wk2i *	x0r	- wk2r * x0i;
			a[j2 + 1] =	-wk2i *	x0i	+ wk2r * x0r;
			x0r	= x1r -	x3i;
			x0i	= x1i +	x3r;
			a[j1] =	wk1r * x0r - wk1i *	x0i;
			a[j1 + 1] =	wk1r * x0i + wk1i *	x0r;
			x0r	= x1r +	x3i;
			x0i	= x1i -	x3r;
			a[j3] =	wk3r * x0r - wk3i *	x0i;
			a[j3 + 1] =	wk3r * x0i + wk3i *	x0r;
		}
	}
}

void rftfsub(int n, float *a, int nc, const float *c)
{
    int j, k, kk, ks, m;
    float wkr, wki, xr, xi, yr, yi;

    m = n >> 1;
    ks = 2 * nc / m;
    kk = 0;
    for (j = 2; j < m; j += 2) {
        k = n - j;
        kk += ks;
        wkr = 0.5 - c[nc - kk];
        wki = c[kk];
        xr = a[j] - a[k];
        xi = a[j + 1] + a[k + 1];
        yr = wkr * xr - wki * xi;
        yi = wkr * xi + wki * xr;
        a[j] -= yr;
        a[j + 1] -= yi;
        a[k] += yr;
        a[k + 1] -= yi;
    }
}

void rftbsub(int n, float *a, int nc, const float *c)
{
    int j, k, kk, ks, m;
    float wkr, wki, xr, xi, yr, yi;

    a[1] = -a[1];
    m = n >> 1;
    ks = 2 * nc / m;
    kk = 0;
    for (j = 2; j < m; j += 2) {
        k = n - j;
        kk += ks;
        wkr = 0.5 - c[nc - kk];
        wki = c[kk];
        xr = a[j] - a[k];
        xi = a[j + 1] + a[k + 1];
        yr = wkr * xr + wki * xi;
        yi = wkr * xi - wki * xr;
        a[j] -= yr;
        a[j + 1] = yi - a[j + 1];
        a[k] += yr;
        a[k + 1] = yi - a[k + 1];
    }
    a[m + 1] = -a[m + 1];
}
