PIFFT_float

git-svn-id: svn://db.shs.com.ru/pip@124 12ceb7fc-bf1f-11e4-8940-5bc7170c53b5
2015-05-04 09:54:02 +00:00
parent f51b36dfed
commit f0d4d75f65
2 changed files with 995 additions and 32 deletions
--- a/src/math/pifft.cpp
+++ b/src/math/pifft.cpp
@@ -20,12 +20,11 @@
 #include "pifft.h"


-PIFFT::PIFFT() {
-	prepared = false;
+PIFFT_double::PIFFT_double() {
 }


-PIVector<complexd> * PIFFT::calcFFT(const PIVector<complexd> & val) {
+PIVector<complexd> * PIFFT_double::calcFFT(const PIVector<complexd> & val) {
 	result.clear();
 	if (val.size_s() < 4) return &result;
 	fftc1d(val, val.size());
@@ -33,7 +32,7 @@ PIVector<complexd> * PIFFT::calcFFT(const PIVector<complexd> & val) {
 }


-PIVector<complexd> * PIFFT::calcFFTinverse(const PIVector<complexd> & val) {
+PIVector<complexd> * PIFFT_double::calcFFTinverse(const PIVector<complexd> & val) {
 	result.clear();
 	if (val.size_s() < 4) return &result;
 	fftc1dinv(val, val.size());
@@ -41,7 +40,7 @@ PIVector<complexd> * PIFFT::calcFFTinverse(const PIVector<complexd> & val) {
 }


-PIVector<complexd> * PIFFT::calcHilbert(const PIVector<double> & val) {
+PIVector<complexd> * PIFFT_double::calcHilbert(const PIVector<double> & val) {
 	result.clear();
 	if (val.size_s() < 4) return &result;
 	fftc1r(val, val.size());
@@ -52,7 +51,7 @@ PIVector<complexd> * PIFFT::calcHilbert(const PIVector<double> & val) {
 }


-PIVector<complexd> * PIFFT::calcFFT(const PIVector<double> & val) {
+PIVector<complexd> * PIFFT_double::calcFFT(const PIVector<double> & val) {
 	result.clear();
 	if (val.size_s() < 4) return &result;
 	fftc1r(val, val.size());
@@ -60,7 +59,7 @@ PIVector<complexd> * PIFFT::calcFFT(const PIVector<double> & val) {
 }


-PIVector<double> PIFFT::getAmplitude() {
+PIVector<double> PIFFT_double::getAmplitude() {
 	PIVector<double> a;
 	double tmp;
 	for (uint i = 0; i < result.size(); i++) {
@@ -71,7 +70,7 @@ PIVector<double> PIFFT::getAmplitude() {
 }


-void PIFFT::fftc1d(const PIVector<complexd> & a, uint n) {
+void PIFFT_double::fftc1d(const PIVector<complexd> & a, uint n) {
 	createPlan(n);
 	uint i;
 	PIVector<double> buf;
@@ -87,7 +86,7 @@ void PIFFT::fftc1d(const PIVector<complexd> & a, uint n) {
 }


-void PIFFT::fftc1r(const PIVector<double> & a, uint n) {
+void PIFFT_double::fftc1r(const PIVector<double> & a, uint n) {
 	uint i;
 	if (n % 2 == 0) {
 		PIVector<double> buf;
@@ -119,7 +118,7 @@ void PIFFT::fftc1r(const PIVector<double> & a, uint n) {
 }


-void PIFFT::fftc1dinv(const PIVector<complexd> & a, uint n) {
+void PIFFT_double::fftc1dinv(const PIVector<complexd> & a, uint n) {
 	PIVector<complexd> cbuf;
 	cbuf.resize(n);
 	uint i;
@@ -133,18 +132,17 @@ void PIFFT::fftc1dinv(const PIVector<complexd> & a, uint n) {
 }


-void PIFFT::createPlan(uint n) {
+void PIFFT_double::createPlan(uint n) {
 	curplan.plan.clear();
 	curplan.precomputed.clear();
 	curplan.stackbuf.clear();
 	curplan.tmpbuf.clear();
 	if (n < 2) return;
 	ftbasegeneratecomplexfftplan(n, &curplan);
-	prepared = true;
 }


-void PIFFT::ftbasegeneratecomplexfftplan(uint n, ftplan * plan) {
+void PIFFT_double::ftbasegeneratecomplexfftplan(uint n, ftplan * plan) {
 	int planarraysize;
 	int plansize;
 	int precomputedsize;
@@ -196,7 +194,7 @@ PARAMETERS:
 -- ALGLIB --
 	Copyright 01.05.2009 by Bochkanov Sergey
 *************************************************************************/
-void PIFFT::ftbase_ftbasegenerateplanrec(
+void PIFFT_double::ftbase_ftbasegenerateplanrec(
 int n,
 int tasktype,
 ftplan * plan,
@@ -294,7 +292,7 @@ Recurrent subroutine for precomputing FFT plans
 -- ALGLIB --
 	Copyright 01.05.2009 by Bochkanov Sergey
 *************************************************************************/
-void PIFFT::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
+void PIFFT_double::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
 int entryoffset,
 ae_int_t stackptr) {
 	int n1, n2, n, m, offs;
@@ -356,7 +354,7 @@ ae_int_t stackptr) {
 }


-void PIFFT::ftbasefactorize(int n, int * n1, int * n2) {
+void PIFFT_double::ftbasefactorize(int n, int * n1, int * n2) {
 	*n1 = *n2 = 0;
 	int ftbase_ftbasecodeletrecommended = 5;
 	if ((*n1) * (*n2) != n) {
@@ -394,7 +392,7 @@ Is number smooth?
 -- ALGLIB --
 	Copyright 01.05.2009 by Bochkanov Sergey
 *************************************************************************/
-void PIFFT::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
+void PIFFT_double::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
 	if (seed >= n) {
 		*best = piMini(*best, seed);
 		return;
@@ -408,7 +406,7 @@ void PIFFT::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * b
 }


-int PIFFT::ftbasefindsmooth(int n) {
+int PIFFT_double::ftbasefindsmooth(int n) {
 	int best, result;
 	best = 2;
 	while (best < n)
@@ -419,13 +417,13 @@ int PIFFT::ftbasefindsmooth(int n) {
 }


-void PIFFT::ftbase_internalreallintranspose(PIVector<double> * a, int m, int n, int astart, PIVector<double> * buf) {
+void PIFFT_double::ftbase_internalreallintranspose(PIVector<double> * a, int m, int n, int astart, PIVector<double> * buf) {
 	ftbase_fftirltrec(a, astart, n, buf, 0, m, m, n);
 	for (int i = 0; i < 2 * m * n; i++) (*a)[astart + i] = (*buf)[i];
 }


-void PIFFT::ftbase_fftirltrec(PIVector<double> * a, int astart, int astride, PIVector<double> * b, int bstart, int bstride, int m, int n) {
+void PIFFT_double::ftbase_fftirltrec(PIVector<double> * a, int astart, int astride, PIVector<double> * b, int bstart, int bstride, int m, int n) {
 	int idx1, idx2;
 	int m1, n1;
 	if (m == 0 || n == 0)
@@ -458,14 +456,14 @@ void PIFFT::ftbase_fftirltrec(PIVector<double> * a, int astart, int astride, PIV
 }


-void PIFFT::ftbase_internalcomplexlintranspose(PIVector<double> * a, int m, int n, int astart, PIVector<double> * buf) {
+void PIFFT_double::ftbase_internalcomplexlintranspose(PIVector<double> * a, int m, int n, int astart, PIVector<double> * buf) {
 	ftbase_ffticltrec(a, astart, n, buf, 0, m, m, n);
 	for (int i = 0; i < 2 * m * n; i++)
 		(*a)[astart + i] = (*buf)[i];
 }


-void PIFFT::ftbase_ffticltrec(PIVector<double> * a, int astart, int astride, PIVector<double> * b, int bstart, int bstride, int m, int n) {
+void PIFFT_double::ftbase_ffticltrec(PIVector<double> * a, int astart, int astride, PIVector<double> * b, int bstart, int bstride, int m, int n) {
 	int idx1, idx2, m2, m1, n1;
 	if (m == 0 || n == 0)
 		return;
@@ -499,7 +497,7 @@ void PIFFT::ftbase_ffticltrec(PIVector<double> * a, int astart, int astride, PIV
 }


-void PIFFT::ftbaseexecuteplan(PIVector<double> * a, int aoffset, int n, ftplan * plan) {
+void PIFFT_double::ftbaseexecuteplan(PIVector<double> * a, int aoffset, int n, ftplan * plan) {
 	ae_int_t stackptr;
 	stackptr = 0;
 	ftbaseexecuteplanrec(a, aoffset, plan, 0, stackptr);
@@ -516,7 +514,7 @@ Parameters:
 -- ALGLIB --
 	Copyright 01.05.2009 by Bochkanov Sergey
 *************************************************************************/
-void PIFFT::ftbaseexecuteplanrec(PIVector<double> * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) {
+void PIFFT_double::ftbaseexecuteplanrec(PIVector<double> * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) {
 	int n1, n2, n, m, offs, offs1, offs2, offsa, offsb, offsp;
 	double hk, hnk, x, y, bx, by, v0, v1, v2, v3;
 	double a0x, a0y, a1x, a1y, a2x, a2y, a3x, a3y;
@@ -884,7 +882,7 @@ Twiddle factors calculation
 -- ALGLIB --
 	Copyright 01.05.2009 by Bochkanov Sergey
 *************************************************************************/
-void PIFFT::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2) {
+void PIFFT_double::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2) {
 	int n, idx, offs;
 	double x, y, twxm1, twy, twbasexm1, twbasey, twrowxm1, twrowy, tmpx, tmpy, v;
 	int ftbase_ftbaseupdatetw = 4;
@@ -934,3 +932,925 @@ void PIFFT::ftbase_ffttwcalc(PIVector<double> * a, int aoffset, int n1, int n2)
 		}
 	}
 }
+
+
+//=================================================================================================
+//*************************************************************************************************
+//---------------------FLOAT-----------------------FLOAT-------------------------------------------
+//*************************************************************************************************
+//=================================================================================================
+
+
+PIFFT_float::PIFFT_float() {
+}
+
+
+PIVector<complexf> * PIFFT_float::calcFFT(const PIVector<complexf> & val) {
+	result.clear();
+	if (val.size_s() < 4) return &result;
+	fftc1d(val, val.size());
+	return &result;
+}
+
+
+PIVector<complexf> * PIFFT_float::calcFFTinverse(const PIVector<complexf> & val) {
+	result.clear();
+	if (val.size_s() < 4) return &result;
+	fftc1dinv(val, val.size());
+	return &result;
+}
+
+
+PIVector<complexf> * PIFFT_float::calcHilbert(const PIVector<float> & val) {
+	result.clear();
+	if (val.size_s() < 4) return &result;
+	fftc1r(val, val.size());
+	for (uint i = 0; i < result.size() / 2; i++)             result[i] = result[i] * 2.f;
+	for (uint i = result.size() / 2; i < result.size(); i++) result[i] = 0;
+	fftc1dinv(result, result.size());
+	return &result;
+}
+
+
+PIVector<complexf> * PIFFT_float::calcFFT(const PIVector<float> & val) {
+	result.clear();
+	if (val.size_s() < 4) return &result;
+	fftc1r(val, val.size());
+	return &result;
+}
+
+
+PIVector<float> PIFFT_float::getAmplitude() {
+	PIVector<float> a;
+	float tmp;
+	for (uint i = 0; i < result.size(); i++) {
+		tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
+		a.push_back(tmp);
+	}
+	return a;
+}
+
+
+void PIFFT_float::fftc1d(const PIVector<complexf> & a, uint n) {
+	createPlan(n);
+	uint i;
+	PIVector<float> buf;
+	buf.resize(2 * n);
+	for (i = 0; i < n; i++) {
+		buf[2 * i + 0] = a[i].real();
+		buf[2 * i + 1] = a[i].imag();
+	}
+	ftbaseexecuteplan(&buf, 0, n, &curplan);
+	result.resize(n);
+	for (i = 0; i < n; i++)
+		result[i] = complexf(buf[2 * i + 0], buf[2 * i + 1]);
+}
+
+
+void PIFFT_float::fftc1r(const PIVector<float> & a, uint n) {
+	uint i;
+	if (n % 2 == 0) {
+		PIVector<float> buf;
+		uint n2 = n / 2;
+		buf = a;
+		createPlan(n2);
+		ftbaseexecuteplan(&buf, 0, n2, &curplan);
+		result.resize(n);
+		uint idx;
+		complexf hn, hmnc, v;
+		for (i = 0; i <= n2; i++) {
+			idx = 2 * (i % n2);
+			hn = complexf(buf[idx + 0], buf[idx + 1]);
+			idx = 2 * ((n2 - i) % n2);
+			hmnc = complexf(buf[idx + 0], -buf[idx + 1]);
+			v = complexf(sin(M_PI * i / n2), cos(M_PI * i / n2));
+			result[i] = ((hn + hmnc) - (v * (hn - hmnc)));
+			result[i] *= 0.5;
+		}
+		for (i = n2 + 1; i < n; i++)
+			result[i] = conj(result[n - i]);
+	} else {
+		PIVector<complexf> cbuf;
+		cbuf.resize(n);
+		for (i = 0; i < n; i++)
+			cbuf[i] = complexf(a[i], 0.);
+		fftc1d(cbuf, n);
+	}
+}
+
+
+void PIFFT_float::fftc1dinv(const PIVector<complexf> & a, uint n) {
+	PIVector<complexf> cbuf;
+	cbuf.resize(n);
+	uint i;
+	for (i = 0; i < n; i++) {
+		cbuf[i] = conj(a[i]);
+	}
+	fftc1d(cbuf, n);
+	for (i = 0; i < n; i++) {
+		result[i] = conj(result[i] / (float)n);
+	}
+}
+
+
+void PIFFT_float::createPlan(uint n) {
+	curplan.plan.clear();
+	curplan.precomputed.clear();
+	curplan.stackbuf.clear();
+	curplan.tmpbuf.clear();
+	if (n < 2) return;
+	ftbasegeneratecomplexfftplan(n, &curplan);
+}
+
+
+void PIFFT_float::ftbasegeneratecomplexfftplan(uint n, ftplan * plan) {
+	int planarraysize;
+	int plansize;
+	int precomputedsize;
+	int tmpmemsize;
+	int stackmemsize;
+	ae_int_t stackptr;
+	planarraysize = 1;
+	plansize = 0;
+	precomputedsize = 0;
+	stackmemsize = 0;
+	stackptr = 0;
+	tmpmemsize = 2 * n;
+	curplan.plan.resize(planarraysize);
+	int ftbase_ftbasecffttask = 0;
+	ftbase_ftbasegenerateplanrec(n, ftbase_ftbasecffttask, plan, &plansize, &precomputedsize, &planarraysize, &tmpmemsize, &stackmemsize, stackptr);
+	if (stackptr != 0) {
+		return;//ae_assert(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
+	}
+	curplan.stackbuf.resize(piMax<int>(stackmemsize, 1)); //ae_vector_set_length(&curplan.stackbuf, ae_maxint(stackmemsize, 1));
+	curplan.tmpbuf.resize(piMax<int>(tmpmemsize, 1)); //ae_vector_set_length(&(curplan.tmpbuf), ae_maxint(tmpmemsize, 1));
+	curplan.precomputed.resize(piMax<int>(precomputedsize, 1)); //ae_vector_set_length(&curplan.precomputed, ae_maxint(precomputedsize, 1));
+	stackptr = 0;
+	ftbase_ftbaseprecomputeplanrec(plan, 0, stackptr);
+	if (stackptr != 0) {
+		return;//ae_assert(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
+	}
+}
+
+
+/*************************************************************************
+Recurrent subroutine for the FFTGeneratePlan:
+
+PARAMETERS:
+	N                   plan size
+	IsReal              whether input is real or not.
+						subroutine MUST NOT ignore this flag because real
+						inputs comes with non-initialized imaginary parts,
+						so ignoring this flag will result in corrupted output
+	HalfOut             whether full output or only half of it from 0 to
+						floor(N/2) is needed. This flag may be ignored if
+						doing so will simplify calculations
+	Plan                plan array
+	PlanSize            size of used part (in integers)
+	PrecomputedSize     size of precomputed array allocated yet
+	PlanArraySize       plan array size (actual)
+	TmpMemSize          temporary memory required size
+	BluesteinMemSize    temporary memory required size
+
+-- ALGLIB --
+	Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void PIFFT_float::ftbase_ftbasegenerateplanrec(
+int n,
+int tasktype,
+ftplan * plan,
+int * plansize,
+int * precomputedsize,
+int * planarraysize,
+int * tmpmemsize,
+int * stackmemsize,
+ae_int_t stackptr, int debugi) {
+	int k, m, n1, n2, esize, entryoffset;
+	int ftbase_ftbaseplanentrysize = 4;
+	int ftbase_ftbasecffttask = 0;
+	int ftbase_fftcooleytukeyplan = 0;
+	int ftbase_fftbluesteinplan = 1;
+	int ftbase_fftcodeletplan = 2;
+	int ftbase_fftrealcooleytukeyplan = 5;
+	int ftbase_fftemptyplan = 6;
+	if (*plansize + ftbase_ftbaseplanentrysize > (*planarraysize)) {
+		curplan.plan.resize(4 * (*planarraysize));
+		*planarraysize = 4 * (*planarraysize);
+	}
+	entryoffset = *plansize;
+	esize = ftbase_ftbaseplanentrysize;
+	*plansize = *plansize + esize;
+	if (n == 1) {
+		curplan.plan[entryoffset + 0] = esize;
+		curplan.plan[entryoffset + 1] = -1;
+		curplan.plan[entryoffset + 2] = -1;
+		curplan.plan[entryoffset + 3] = ftbase_fftemptyplan;
+		curplan.plan[entryoffset + 4] = -1;
+		curplan.plan[entryoffset + 5] = -1;
+		curplan.plan[entryoffset + 6] = -1;
+		curplan.plan[entryoffset + 7] = -1;
+		return;
+	}
+	ftbasefactorize(n, &n1, &n2);
+	if (n1 != 1) {
+		*tmpmemsize = piMax<int>(*tmpmemsize, 2 * n1 * n2);
+		curplan.plan[entryoffset + 0] = esize;
+		curplan.plan[entryoffset + 1] = n1;
+		curplan.plan[entryoffset + 2] = n2;
+		if (tasktype == ftbase_ftbasecffttask)
+			curplan.plan[entryoffset + 3] = ftbase_fftcooleytukeyplan;
+		else
+			curplan.plan[entryoffset + 3] = ftbase_fftrealcooleytukeyplan;
+		curplan.plan[entryoffset + 4] = 0;
+		curplan.plan[entryoffset + 5] = *plansize;
+		debugi++;
+		ftbase_ftbasegenerateplanrec(n1, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr, debugi);
+		curplan.plan[entryoffset + 6] = *plansize;
+		ftbase_ftbasegenerateplanrec(n2, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr, debugi);
+		curplan.plan[entryoffset + 7] = -1;
+		return;
+	} else {
+		if (n >= 2 && n <= 5) {
+			curplan.plan[entryoffset + 0] = esize;
+			curplan.plan[entryoffset + 1] = n1;
+			curplan.plan[entryoffset + 2] = n2;
+			curplan.plan[entryoffset + 3] = ftbase_fftcodeletplan;
+			curplan.plan[entryoffset + 4] = 0;
+			curplan.plan[entryoffset + 5] = -1;
+			curplan.plan[entryoffset + 6] = -1;
+			curplan.plan[entryoffset + 7] = *precomputedsize;
+			if (n == 3)
+				*precomputedsize = *precomputedsize + 2;
+			if (n == 5)
+				*precomputedsize = *precomputedsize + 5;
+			return;
+		} else {
+			k = 2 * n2 - 1;
+			m = ftbasefindsmooth(k);
+			*tmpmemsize = piMax<int>(*tmpmemsize, 2 * m);
+			curplan.plan[entryoffset + 0] = esize;
+			curplan.plan[entryoffset + 1] = n2;
+			curplan.plan[entryoffset + 2] = -1;
+			curplan.plan[entryoffset + 3] = ftbase_fftbluesteinplan;
+			curplan.plan[entryoffset + 4] = m;
+			curplan.plan[entryoffset + 5] = *plansize;
+			stackptr = stackptr + 2 * 2 * m;
+			*stackmemsize = piMax<int>(*stackmemsize, stackptr);
+			ftbase_ftbasegenerateplanrec(m, ftbase_ftbasecffttask, plan, plansize, precomputedsize, planarraysize, tmpmemsize, stackmemsize, stackptr);
+			stackptr = stackptr - 2 * 2 * m;
+			curplan.plan[entryoffset + 6] = -1;
+			curplan.plan[entryoffset + 7] = *precomputedsize;
+			*precomputedsize = *precomputedsize + 2 * m + 2 * n;
+			return;
+		}
+	}
+}
+
+
+/*************************************************************************
+Recurrent subroutine for precomputing FFT plans
+
+-- ALGLIB --
+	Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void PIFFT_float::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
+int entryoffset,
+ae_int_t stackptr) {
+	int n1, n2, n, m, offs;
+	float v, bx, by;
+	int ftbase_fftcooleytukeyplan = 0;
+	int ftbase_fftbluesteinplan = 1;
+	int ftbase_fftcodeletplan = 2;
+	int ftbase_fhtcooleytukeyplan = 3;
+	int ftbase_fhtcodeletplan = 4;
+	int ftbase_fftrealcooleytukeyplan = 5;
+	if ((curplan.plan[entryoffset + 3] == ftbase_fftcooleytukeyplan || curplan.plan[entryoffset + 3] == ftbase_fftrealcooleytukeyplan) || curplan.plan[entryoffset + 3] == ftbase_fhtcooleytukeyplan) {
+		ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 5], stackptr);
+		ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 6], stackptr);
+		return;
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fftcodeletplan || curplan.plan[entryoffset + 3] == ftbase_fhtcodeletplan) {
+		n1 = curplan.plan[entryoffset + 1];
+		n2 = curplan.plan[entryoffset + 2];
+		n = n1 * n2;
+		if (n == 3) {
+			offs = curplan.plan[entryoffset + 7];
+			curplan.precomputed[offs + 0] = cos(2 * M_PI / 3) - 1;
+			curplan.precomputed[offs + 1] = sin(2 * M_PI / 3);
+			return;
+		}
+		if (n == 5) {
+			offs = curplan.plan[entryoffset + 7];
+			v = 2 * M_PI / 5;
+			curplan.precomputed[offs + 0] = (cos(v) + cos(2 * v)) / 2 - 1;
+			curplan.precomputed[offs + 1] = (cos(v) - cos(2 * v)) / 2;
+			curplan.precomputed[offs + 2] = -sin(v);
+			curplan.precomputed[offs + 3] = -(sin(v) + sin(2 * v));
+			curplan.precomputed[offs + 4] = sin(v) - sin(2 * v);
+			return;
+		}
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fftbluesteinplan) {
+		ftbase_ftbaseprecomputeplanrec(plan, curplan.plan[entryoffset + 5], stackptr);
+		n = curplan.plan[entryoffset + 1];
+		m = curplan.plan[entryoffset + 4];
+		offs = curplan.plan[entryoffset + 7];
+		for (int i = 0; i <= 2 * m - 1; i++)
+			curplan.precomputed[offs + i] = 0;
+		for (int i = 0; i < n; i++) {
+			bx = cos(M_PI * sqr(i) / n);
+			by = sin(M_PI * sqr(i) / n);
+			curplan.precomputed[offs + 2 * i + 0] = bx;
+			curplan.precomputed[offs + 2 * i + 1] = by;
+			curplan.precomputed[offs + 2 * m + 2 * i + 0] = bx;
+			curplan.precomputed[offs + 2 * m + 2 * i + 1] = by;
+			if (i > 0) {
+				curplan.precomputed[offs + 2 * (m - i) + 0] = bx;
+				curplan.precomputed[offs + 2 * (m - i) + 1] = by;
+			}
+		}
+		ftbaseexecuteplanrec(&curplan.precomputed, offs, plan, curplan.plan[entryoffset + 5], stackptr);
+		return;
+	}
+}
+
+
+void PIFFT_float::ftbasefactorize(int n, int * n1, int * n2) {
+	*n1 = *n2 = 0;
+	int ftbase_ftbasecodeletrecommended = 5;
+	if ((*n1) * (*n2) != n) {
+		for (int j = ftbase_ftbasecodeletrecommended; j >= 2; j--) {
+			if (n % j == 0) {
+				*n1 = j;
+				*n2 = n / j;
+				break;
+			}
+		}
+	}
+	if ((*n1) * (*n2) != n) {
+		for (int j = ftbase_ftbasecodeletrecommended + 1; j <= n - 1; j++) {
+			if (n % j == 0) {
+				*n1 = j;
+				*n2 = n / j;
+				break;
+			}
+		}
+	}
+	if ((*n1) * (*n2) != n) {
+		*n1 = 1;
+		*n2 = n;
+	}
+	if ((*n2) == 1 && (*n1) != 1) {
+		*n2 = *n1;
+		*n1 = 1;
+	}
+}
+
+
+/*************************************************************************
+Is number smooth?
+
+-- ALGLIB --
+	Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void PIFFT_float::ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int * best) {
+	if (seed >= n) {
+		*best = piMini(*best, seed);
+		return;
+	}
+	if (leastfactor <= 2)
+		ftbase_ftbasefindsmoothrec(n, seed * 2, 2, best);
+	if (leastfactor <= 3)
+		ftbase_ftbasefindsmoothrec(n, seed * 3, 3, best);
+	if (leastfactor <= 5)
+		ftbase_ftbasefindsmoothrec(n, seed * 5, 5, best);
+}
+
+
+int PIFFT_float::ftbasefindsmooth(int n) {
+	int best, result;
+	best = 2;
+	while (best < n)
+		best = 2 * best;
+	ftbase_ftbasefindsmoothrec(n, 1, 2, &best);
+	result = best;
+	return result;
+}
+
+
+void PIFFT_float::ftbase_internalreallintranspose(PIVector<float> * a, int m, int n, int astart, PIVector<float> * buf) {
+	ftbase_fftirltrec(a, astart, n, buf, 0, m, m, n);
+	for (int i = 0; i < 2 * m * n; i++) (*a)[astart + i] = (*buf)[i];
+}
+
+
+void PIFFT_float::ftbase_fftirltrec(PIVector<float> * a, int astart, int astride, PIVector<float> * b, int bstart, int bstride, int m, int n) {
+	int idx1, idx2;
+	int m1, n1;
+	if (m == 0 || n == 0)
+		return;
+	if (piMaxi(m, n) <= 8) {
+		for (int i = 0; i <= m - 1; i++) {
+			idx1 = bstart + i;
+			idx2 = astart + i * astride;
+			for (int j = 0; j <= n - 1; j++) {
+				(*b)[idx1] = a->at(idx2);
+				idx1 = idx1 + bstride;
+				idx2 = idx2 + 1;
+			}
+		}
+		return;
+	}
+	if (n > m) {
+		n1 = n / 2;
+		if (n - n1 >= 8 && n1 % 8 != 0)
+			n1 = n1 + (8 - n1 % 8);
+		ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m, n1);
+		ftbase_fftirltrec(a, astart + n1, astride, b, bstart + n1 * bstride, bstride, m, n - n1);
+	} else {
+		m1 = m / 2;
+		if (m - m1 >= 8 && m1 % 8 != 0)
+			m1 = m1 + (8 - m1 % 8);
+		ftbase_fftirltrec(a, astart, astride, b, bstart, bstride, m1, n);
+		ftbase_fftirltrec(a, astart + m1 * astride, astride, b, bstart + m1, bstride, m - m1, n);
+	}
+}
+
+
+void PIFFT_float::ftbase_internalcomplexlintranspose(PIVector<float> * a, int m, int n, int astart, PIVector<float> * buf) {
+	ftbase_ffticltrec(a, astart, n, buf, 0, m, m, n);
+	for (int i = 0; i < 2 * m * n; i++)
+		(*a)[astart + i] = (*buf)[i];
+}
+
+
+void PIFFT_float::ftbase_ffticltrec(PIVector<float> * a, int astart, int astride, PIVector<float> * b, int bstart, int bstride, int m, int n) {
+	int idx1, idx2, m2, m1, n1;
+	if (m == 0 || n == 0)
+		return;
+	if (piMax<int>(m, n) <= 8) {
+		m2 = 2 * bstride;
+		for (int i = 0; i <= m - 1; i++) {
+			idx1 = bstart + 2 * i;
+			idx2 = astart + 2 * i * astride;
+			for (int j = 0; j <= n - 1; j++) {
+				(*b)[idx1 + 0] = a->at(idx2 + 0);
+				(*b)[idx1 + 1] = a->at(idx2 + 1);
+				idx1 = idx1 + m2;
+				idx2 = idx2 + 2;
+			}
+		}
+		return;
+	}
+	if (n > m) {
+		n1 = n / 2;
+		if (n - n1 >= 8 && n1 % 8 != 0)
+			n1 = n1 + (8 - n1 % 8);
+		ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m, n1);
+		ftbase_ffticltrec(a, astart + 2 * n1, astride, b, bstart + 2 * n1 * bstride, bstride, m, n - n1);
+	} else {
+		m1 = m / 2;
+		if (m - m1 >= 8 && m1 % 8 != 0)
+			m1 = m1 + (8 - m1 % 8);
+		ftbase_ffticltrec(a, astart, astride, b, bstart, bstride, m1, n);
+		ftbase_ffticltrec(a, astart + 2 * m1 * astride, astride, b, bstart + 2 * m1, bstride, m - m1, n);
+	}
+}
+
+
+void PIFFT_float::ftbaseexecuteplan(PIVector<float> * a, int aoffset, int n, ftplan * plan) {
+	ae_int_t stackptr;
+	stackptr = 0;
+	ftbaseexecuteplanrec(a, aoffset, plan, 0, stackptr);
+}
+
+
+/*************************************************************************
+Recurrent subroutine for the FTBaseExecutePlan
+
+Parameters:
+	A           FFT'ed array
+	AOffset     offset of the FFT'ed part (distance is measured in floats)
+
+-- ALGLIB --
+	Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void PIFFT_float::ftbaseexecuteplanrec(PIVector<float> * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) {
+	int n1, n2, n, m, offs, offs1, offs2, offsa, offsb, offsp;
+	float hk, hnk, x, y, bx, by, v0, v1, v2, v3;
+	float a0x, a0y, a1x, a1y, a2x, a2y, a3x, a3y;
+	float t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, t5x, t5y;
+	float m1x, m1y, m2x, m2y, m3x, m3y, m4x, m4y, m5x, m5y;
+	float s1x, s1y, s2x, s2y, s3x, s3y, s4x, s4y, s5x, s5y;
+	float c1, c2, c3, c4, c5;
+	int ftbase_fftcooleytukeyplan = 0;
+	int ftbase_fftbluesteinplan = 1;
+	int ftbase_fftcodeletplan = 2;
+	int ftbase_fhtcooleytukeyplan = 3;
+	int ftbase_fhtcodeletplan = 4;
+	int ftbase_fftrealcooleytukeyplan = 5;
+	int ftbase_fftemptyplan = 6;
+	PIVector<float> & tmpb(curplan.tmpbuf);
+
+	if (curplan.plan[entryoffset + 3] == ftbase_fftemptyplan)
+		return;
+	if (curplan.plan[entryoffset + 3] == ftbase_fftcooleytukeyplan) {
+		n1 = curplan.plan[entryoffset + 1];
+		n2 = curplan.plan[entryoffset + 2];
+		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
+		for (int i = 0; i <= n2 - 1; i++)
+			ftbaseexecuteplanrec(a, aoffset + i * n1 * 2, plan, curplan.plan[entryoffset + 5], stackptr);
+		ftbase_ffttwcalc(a, aoffset, n1, n2);
+		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
+		for (int i = 0; i <= n1 - 1; i++)
+			ftbaseexecuteplanrec(a, aoffset + i * n2 * 2, plan, curplan.plan[entryoffset + 6], stackptr);
+		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
+		return;
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fftrealcooleytukeyplan) {
+		n1 = curplan.plan[entryoffset + 1];
+		n2 = curplan.plan[entryoffset + 2];
+		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
+		for (int i = 0; i <= n1 / 2 - 1; i++) {
+			offs = aoffset + 2 * i * n2 * 2;
+			for (int k = 0; k <= n2 - 1; k++)
+				(*a)[offs + 2 * k + 1] = (*a)[offs + 2 * n2 + 2 * k + 0];
+			ftbaseexecuteplanrec(a, offs, plan, curplan.plan[entryoffset + 6], stackptr);
+			tmpb[0] = (*a)[offs + 0];
+			tmpb[1] = 0;
+			tmpb[2 * n2 + 0] = (*a)[offs + 1];
+			tmpb[2 * n2 + 1] = 0;
+			for (int k = 1; k <= n2 - 1; k++) {
+				offs1 = 2 * k;
+				offs2 = 2 * n2 + 2 * k;
+				hk = (*a)[offs + 2 * k + 0];
+				hnk = (*a)[offs + 2 * (n2 - k) + 0];
+				tmpb[offs1 + 0] = 0.5 * (hk + hnk);
+				tmpb[offs2 + 1] = -0.5 * (hk - hnk);
+				hk = (*a)[offs + 2 * k + 1];
+				hnk = (*a)[offs + 2 * (n2 - k) + 1];
+				tmpb[offs2 + 0] = 0.5 * (hk + hnk);
+				tmpb[offs1 + 1] = 0.5 * (hk - hnk);
+			}
+			for (int k = 0; k < 2 * n2 * 2; k++) (*a)[offs + k] = tmpb[k];
+		}
+		if (n1 % 2 != 0)
+			ftbaseexecuteplanrec(a, aoffset + (n1 - 1)*n2 * 2, plan, curplan.plan[entryoffset + 6], stackptr);
+		ftbase_ffttwcalc(a, aoffset, n2, n1);
+		ftbase_internalcomplexlintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
+		for (int i = 0; i <= n2 - 1; i++)
+			ftbaseexecuteplanrec(a, aoffset + i * n1 * 2, plan, curplan.plan[entryoffset + 5], stackptr);
+		ftbase_internalcomplexlintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
+		return;
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fhtcooleytukeyplan) {
+		n1 = curplan.plan[entryoffset + 1];
+		n2 = curplan.plan[entryoffset + 2];
+		n = n1 * n2;
+		ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
+		for (int i = 0; i <= n2 - 1; i++)
+			ftbaseexecuteplanrec(a, aoffset + i * n1, plan, curplan.plan[entryoffset + 5], stackptr);
+		for (int i = 0; i <= n2 - 1; i++) {
+			for (int j = 0; j <= n1 - 1; j++) {
+				offsa = aoffset + i * n1;
+				hk = (*a)[offsa + j];
+				hnk = (*a)[offsa + (n1 - j) % n1];
+				offs = 2 * (i * n1 + j);
+				tmpb[offs + 0] = -0.5 * (hnk - hk);
+				tmpb[offs + 1] = 0.5 * (hk + hnk);
+			}
+		}
+		ftbase_ffttwcalc(&(curplan.tmpbuf), 0, n1, n2);
+		for (int j = 0; j <= n1 - 1; j++)
+			(*a)[aoffset + j] = tmpb[2 * j + 0] + tmpb[2 * j + 1];
+		if (n2 % 2 == 0) {
+			offs = 2 * (n2 / 2) * n1;
+			offsa = aoffset + n2 / 2 * n1;
+			for (int j = 0; j <= n1 - 1; j++)
+				(*a)[offsa + j] = tmpb[offs + 2 * j + 0] + tmpb[offs + 2 * j + 1];
+		}
+		for (int i = 1; i <= (n2 + 1) / 2 - 1; i++) {
+			offs = 2 * i * n1;
+			offs2 = 2 * (n2 - i) * n1;
+			offsa = aoffset + i * n1;
+			for (int j = 0; j <= n1 - 1; j++)
+				(*a)[offsa + j] = tmpb[offs + 2 * j + 1] + tmpb[offs2 + 2 * j + 0];
+			offsa = aoffset + (n2 - i) * n1;
+			for (int j = 0; j <= n1 - 1; j++)
+				(*a)[offsa + j] = tmpb[offs + 2 * j + 0] + tmpb[offs2 + 2 * j + 1];
+		}
+		ftbase_internalreallintranspose(a, n2, n1, aoffset, &(curplan.tmpbuf));
+		for (int i = 0; i <= n1 - 1; i++)
+			ftbaseexecuteplanrec(a, aoffset + i * n2, plan, curplan.plan[entryoffset + 6], stackptr);
+		ftbase_internalreallintranspose(a, n1, n2, aoffset, &(curplan.tmpbuf));
+		return;
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fftcodeletplan) {
+		n1 = curplan.plan[entryoffset + 1];
+		n2 = curplan.plan[entryoffset + 2];
+		n = n1 * n2;
+		if (n == 2) {
+			a0x = (*a)[aoffset + 0];
+			a0y = (*a)[aoffset + 1];
+			a1x = (*a)[aoffset + 2];
+			a1y = (*a)[aoffset + 3];
+			v0 = a0x + a1x;
+			v1 = a0y + a1y;
+			v2 = a0x - a1x;
+			v3 = a0y - a1y;
+			(*a)[aoffset + 0] = v0;
+			(*a)[aoffset + 1] = v1;
+			(*a)[aoffset + 2] = v2;
+			(*a)[aoffset + 3] = v3;
+			return;
+		}
+		if (n == 3) {
+			offs = curplan.plan[entryoffset + 7];
+			c1 = curplan.precomputed[offs + 0];
+			c2 = curplan.precomputed[offs + 1];
+			a0x = (*a)[aoffset + 0];
+			a0y = (*a)[aoffset + 1];
+			a1x = (*a)[aoffset + 2];
+			a1y = (*a)[aoffset + 3];
+			a2x = (*a)[aoffset + 4];
+			a2y = (*a)[aoffset + 5];
+			t1x = a1x + a2x;
+			t1y = a1y + a2y;
+			a0x = a0x + t1x;
+			a0y = a0y + t1y;
+			m1x = c1 * t1x;
+			m1y = c1 * t1y;
+			m2x = c2 * (a1y - a2y);
+			m2y = c2 * (a2x - a1x);
+			s1x = a0x + m1x;
+			s1y = a0y + m1y;
+			a1x = s1x + m2x;
+			a1y = s1y + m2y;
+			a2x = s1x - m2x;
+			a2y = s1y - m2y;
+			(*a)[aoffset + 0] = a0x;
+			(*a)[aoffset + 1] = a0y;
+			(*a)[aoffset + 2] = a1x;
+			(*a)[aoffset + 3] = a1y;
+			(*a)[aoffset + 4] = a2x;
+			(*a)[aoffset + 5] = a2y;
+			return;
+		}
+		if (n == 4) {
+			a0x = (*a)[aoffset + 0];
+			a0y = (*a)[aoffset + 1];
+			a1x = (*a)[aoffset + 2];
+			a1y = (*a)[aoffset + 3];
+			a2x = (*a)[aoffset + 4];
+			a2y = (*a)[aoffset + 5];
+			a3x = (*a)[aoffset + 6];
+			a3y = (*a)[aoffset + 7];
+			t1x = a0x + a2x;
+			t1y = a0y + a2y;
+			t2x = a1x + a3x;
+			t2y = a1y + a3y;
+			m2x = a0x - a2x;
+			m2y = a0y - a2y;
+			m3x = a1y - a3y;
+			m3y = a3x - a1x;
+			(*a)[aoffset + 0] = t1x + t2x;
+			(*a)[aoffset + 1] = t1y + t2y;
+			(*a)[aoffset + 4] = t1x - t2x;
+			(*a)[aoffset + 5] = t1y - t2y;
+			(*a)[aoffset + 2] = m2x + m3x;
+			(*a)[aoffset + 3] = m2y + m3y;
+			(*a)[aoffset + 6] = m2x - m3x;
+			(*a)[aoffset + 7] = m2y - m3y;
+			return;
+		}
+		if (n == 5) {
+			offs = curplan.plan[entryoffset + 7];
+			c1 = curplan.precomputed[offs + 0];
+			c2 = curplan.precomputed[offs + 1];
+			c3 = curplan.precomputed[offs + 2];
+			c4 = curplan.precomputed[offs + 3];
+			c5 = curplan.precomputed[offs + 4];
+			t1x = (*a)[aoffset + 2] + (*a)[aoffset + 8];
+			t1y = (*a)[aoffset + 3] + (*a)[aoffset + 9];
+			t2x = (*a)[aoffset + 4] + (*a)[aoffset + 6];
+			t2y = (*a)[aoffset + 5] + (*a)[aoffset + 7];
+			t3x = (*a)[aoffset + 2] - (*a)[aoffset + 8];
+			t3y = (*a)[aoffset + 3] - (*a)[aoffset + 9];
+			t4x = (*a)[aoffset + 6] - (*a)[aoffset + 4];
+			t4y = (*a)[aoffset + 7] - (*a)[aoffset + 5];
+			t5x = t1x + t2x;
+			t5y = t1y + t2y;
+			(*a)[aoffset + 0] = (*a)[aoffset + 0] + t5x;
+			(*a)[aoffset + 1] = (*a)[aoffset + 1] + t5y;
+			m1x = c1 * t5x;
+			m1y = c1 * t5y;
+			m2x = c2 * (t1x - t2x);
+			m2y = c2 * (t1y - t2y);
+			m3x = -c3 * (t3y + t4y);
+			m3y = c3 * (t3x + t4x);
+			m4x = -c4 * t4y;
+			m4y = c4 * t4x;
+			m5x = -c5 * t3y;
+			m5y = c5 * t3x;
+			s3x = m3x - m4x;
+			s3y = m3y - m4y;
+			s5x = m3x + m5x;
+			s5y = m3y + m5y;
+			s1x = (*a)[aoffset + 0] + m1x;
+			s1y = (*a)[aoffset + 1] + m1y;
+			s2x = s1x + m2x;
+			s2y = s1y + m2y;
+			s4x = s1x - m2x;
+			s4y = s1y - m2y;
+			(*a)[aoffset + 2] = s2x + s3x;
+			(*a)[aoffset + 3] = s2y + s3y;
+			(*a)[aoffset + 4] = s4x + s5x;
+			(*a)[aoffset + 5] = s4y + s5y;
+			(*a)[aoffset + 6] = s4x - s5x;
+			(*a)[aoffset + 7] = s4y - s5y;
+			(*a)[aoffset + 8] = s2x - s3x;
+			(*a)[aoffset + 9] = s2y - s3y;
+			return;
+		}
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fhtcodeletplan) {
+		n1 = curplan.plan[entryoffset + 1];
+		n2 = curplan.plan[entryoffset + 2];
+		n = n1 * n2;
+		if (n == 2) {
+			a0x = (*a)[aoffset + 0];
+			a1x = (*a)[aoffset + 1];
+			(*a)[aoffset + 0] = a0x + a1x;
+			(*a)[aoffset + 1] = a0x - a1x;
+			return;
+		}
+		if (n == 3) {
+			offs = curplan.plan[entryoffset + 7];
+			c1 = curplan.precomputed[offs + 0];
+			c2 = curplan.precomputed[offs + 1];
+			a0x = (*a)[aoffset + 0];
+			a1x = (*a)[aoffset + 1];
+			a2x = (*a)[aoffset + 2];
+			t1x = a1x + a2x;
+			a0x = a0x + t1x;
+			m1x = c1 * t1x;
+			m2y = c2 * (a2x - a1x);
+			s1x = a0x + m1x;
+			(*a)[aoffset + 0] = a0x;
+			(*a)[aoffset + 1] = s1x - m2y;
+			(*a)[aoffset + 2] = s1x + m2y;
+			return;
+		}
+		if (n == 4) {
+			a0x = (*a)[aoffset + 0];
+			a1x = (*a)[aoffset + 1];
+			a2x = (*a)[aoffset + 2];
+			a3x = (*a)[aoffset + 3];
+			t1x = a0x + a2x;
+			t2x = a1x + a3x;
+			m2x = a0x - a2x;
+			m3y = a3x - a1x;
+			(*a)[aoffset + 0] = t1x + t2x;
+			(*a)[aoffset + 1] = m2x - m3y;
+			(*a)[aoffset + 2] = t1x - t2x;
+			(*a)[aoffset + 3] = m2x + m3y;
+			return;
+		}
+		if (n == 5) {
+			offs = curplan.plan[entryoffset + 7];
+			c1 = curplan.precomputed[offs + 0];
+			c2 = curplan.precomputed[offs + 1];
+			c3 = curplan.precomputed[offs + 2];
+			c4 = curplan.precomputed[offs + 3];
+			c5 = curplan.precomputed[offs + 4];
+			t1x = (*a)[aoffset + 1] + (*a)[aoffset + 4];
+			t2x = (*a)[aoffset + 2] + (*a)[aoffset + 3];
+			t3x = (*a)[aoffset + 1] - (*a)[aoffset + 4];
+			t4x = (*a)[aoffset + 3] - (*a)[aoffset + 2];
+			t5x = t1x + t2x;
+			v0 = (*a)[aoffset + 0] + t5x;
+			(*a)[aoffset + 0] = v0;
+			m2x = c2 * (t1x - t2x);
+			m3y = c3 * (t3x + t4x);
+			s3y = m3y - c4 * t4x;
+			s5y = m3y + c5 * t3x;
+			s1x = v0 + c1 * t5x;
+			s2x = s1x + m2x;
+			s4x = s1x - m2x;
+			(*a)[aoffset + 1] = s2x - s3y;
+			(*a)[aoffset + 2] = s4x - s5y;
+			(*a)[aoffset + 3] = s4x + s5y;
+			(*a)[aoffset + 4] = s2x + s3y;
+			return;
+		}
+	}
+	if (curplan.plan[entryoffset + 3] == ftbase_fftbluesteinplan) {
+		n = curplan.plan[entryoffset + 1];
+		m = curplan.plan[entryoffset + 4];
+		offs = curplan.plan[entryoffset + 7];
+		for (int i = stackptr + 2 * n; i <= stackptr + 2 * m - 1; i++)
+			curplan.stackbuf[i] = 0;
+		offsp = offs + 2 * m;
+		offsa = aoffset;
+		offsb = stackptr;
+		for (int i = 0; i < n; i++) {
+			bx = curplan.precomputed[offsp + 0];
+			by = curplan.precomputed[offsp + 1];
+			x = (*a)[offsa + 0];
+			y = (*a)[offsa + 1];
+			curplan.stackbuf[offsb + 0] = x * bx - y * (-by);
+			curplan.stackbuf[offsb + 1] = x * (-by) + y * bx;
+			offsp = offsp + 2;
+			offsa = offsa + 2;
+			offsb = offsb + 2;
+		}
+		ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset + 5], stackptr + 2 * 2 * m);
+		offsb = stackptr;
+		offsp = offs;
+		for (int i = 0; i <= m - 1; i++) {
+			x = curplan.stackbuf[offsb + 0];
+			y = curplan.stackbuf[offsb + 1];
+			bx = curplan.precomputed[offsp + 0];
+			by = curplan.precomputed[offsp + 1];
+			curplan.stackbuf[offsb + 0] = x * bx - y * by;
+			curplan.stackbuf[offsb + 1] = -(x * by + y * bx);
+			offsb = offsb + 2;
+			offsp = offsp + 2;
+		}
+		ftbaseexecuteplanrec(&curplan.stackbuf, stackptr, plan, curplan.plan[entryoffset + 5], stackptr + 2 * 2 * m);
+		offsb = stackptr;
+		offsp = offs + 2 * m;
+		offsa = aoffset;
+		for (int i = 0; i < n; i++) {
+			x = curplan.stackbuf[offsb + 0] / m;
+			y = -curplan.stackbuf[offsb + 1] / m;
+			bx = curplan.precomputed[offsp + 0];
+			by = curplan.precomputed[offsp + 1];
+			(*a)[offsa + 0] = x * bx - y * (-by);
+			(*a)[offsa + 1] = x * (-by) + y * bx;
+			offsp = offsp + 2;
+			offsa = offsa + 2;
+			offsb = offsb + 2;
+		}
+		return;
+	}
+}
+
+
+/*************************************************************************
+Twiddle factors calculation
+
+-- ALGLIB --
+	Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************/
+void PIFFT_float::ftbase_ffttwcalc(PIVector<float> * a, int aoffset, int n1, int n2) {
+	int n, idx, offs;
+	float x, y, twxm1, twy, twbasexm1, twbasey, twrowxm1, twrowy, tmpx, tmpy, v;
+	int ftbase_ftbaseupdatetw = 4;
+	n = n1 * n2;
+	v = -2 * M_PI / n;
+	twbasexm1 = -2 * sqr(sin(0.5 * v));
+	twbasey = sin(v);
+	twrowxm1 = 0;
+	twrowy = 0;
+	for (int i = 0, j = 0; i <= n2 - 1; i++) {
+		twxm1 = 0;
+		twy = 0;
+		for (j = 0; j <= n1 - 1; j++) {
+			idx = i * n1 + j;
+			offs = aoffset + 2 * idx;
+			x = (*a)[offs + 0];
+			y = (*a)[offs + 1];
+			tmpx = x * twxm1 - y * twy;
+			tmpy = x * twy + y * twxm1;
+			(*a)[offs + 0] = x + tmpx;
+			(*a)[offs + 1] = y + tmpy;
+			if (j < n1 - 1) {
+				if (j % ftbase_ftbaseupdatetw == 0) {
+					v = -2 * M_PI * i * (j + 1) / n;
+					twxm1 = -2 * sqr(sin(0.5 * v));
+					twy = sin(v);
+				} else {
+					tmpx = twrowxm1 + twxm1 * twrowxm1 - twy * twrowy;
+					tmpy = twrowy + twxm1 * twrowy + twy * twrowxm1;
+					twxm1 = twxm1 + tmpx;
+					twy = twy + tmpy;
+				}
+			}
+		}
+
+		if (i < n2 - 1) {
+			if (j % ftbase_ftbaseupdatetw == 0) {
+				v = -2 * M_PI * (i + 1) / n;
+				twrowxm1 = -2 * sqr(sin(0.5 * v));
+				twrowy = sin(v);
+			} else {
+				tmpx = twbasexm1 + twrowxm1 * twbasexm1 - twrowy * twbasey;
+				tmpy = twbasey + twrowxm1 * twbasey + twrowy * twbasexm1;
+				twrowxm1 = twrowxm1 + tmpx;
+				twrowy = twrowy + tmpy;
+			}
+		}
+	}
+}
+
--- a/src/math/pifft.h
+++ b/src/math/pifft.h
@@ -25,10 +25,10 @@

 #include "pimathbase.h"

-class PIP_EXPORT PIFFT
+class PIP_EXPORT PIFFT_double
 {
 public:
-	PIFFT();
+	PIFFT_double();
 	
 	PIVector<complexd> * calcFFT(const PIVector<complexd> &val);
 	PIVector<complexd> * calcFFT(const PIVector<double> &val);
@@ -38,11 +38,7 @@ public:
 	
 private:
 	PIVector<complexd> result;
-	bool prepared;
 	typedef ptrdiff_t ae_int_t;
-	void calc_coefs(uint cnt2);
-	void calc_indexes(uint cnt2, uint deep2);
-	complexd coef(uint n, uint k);

 	struct ftplan {
 		PIVector<int> plan;
@@ -71,7 +67,54 @@ private:
 	void ftbase_internalreallintranspose(PIVector<double> *a, int m, int n, int astart, PIVector<double> *buf);
 	void ftbase_fftirltrec(PIVector<double> *a, int astart, int astride, PIVector<double> *b, int bstart, int bstride, int m, int n);
 	void ftbase_ffttwcalc(PIVector<double> *a, int aoffset, int n1, int n2);
-	
 };

+class PIP_EXPORT PIFFT_float
+{
+public:
+	PIFFT_float();
+
+	PIVector<complexf> * calcFFT(const PIVector<complexf> &val);
+	PIVector<complexf> * calcFFT(const PIVector<float> &val);
+	PIVector<complexf> * calcFFTinverse(const PIVector<complexf> &val);
+	PIVector<complexf> * calcHilbert(const PIVector<float> &val);
+	PIVector<float> getAmplitude();
+
+private:
+	PIVector<complexf> result;
+	typedef ptrdiff_t ae_int_t;
+
+	struct ftplan {
+		PIVector<int> plan;
+		PIVector<float> precomputed;
+		PIVector<float> tmpbuf;
+		PIVector<float> stackbuf;
+	};
+
+	ftplan curplan;
+
+	void fftc1d(const PIVector<complexf> &a, uint n);
+	void fftc1r(const PIVector<float> &a, uint n);
+	void fftc1dinv(const PIVector<complexf> &a, uint n);
+
+	void createPlan(uint n);
+	void ftbasegeneratecomplexfftplan(uint n, ftplan *plan);
+	void ftbase_ftbasegenerateplanrec(int n, int tasktype, ftplan *plan, int *plansize, int *precomputedsize, int *planarraysize, int *tmpmemsize, int *stackmemsize, ae_int_t stackptr, int debugi=0);
+	void ftbase_ftbaseprecomputeplanrec(ftplan *plan, int entryoffset, ae_int_t stackptr);
+	void ftbasefactorize(int n, int *n1, int *n2);
+	void ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int *best);
+	int ftbasefindsmooth(int n);
+	void ftbaseexecuteplan(PIVector<float> *a, int aoffset, int n, ftplan *plan);
+	void ftbaseexecuteplanrec(PIVector<float> *a, int aoffset, ftplan *plan, int entryoffset, ae_int_t stackptr);
+	void ftbase_internalcomplexlintranspose(PIVector<float> *a, int m, int n, int astart, PIVector<float> *buf);
+	void ftbase_ffticltrec(PIVector<float> *a, int astart, int astride, PIVector<float> *b, int bstart, int bstride, int m, int n);
+	void ftbase_internalreallintranspose(PIVector<float> *a, int m, int n, int astart, PIVector<float> *buf);
+	void ftbase_fftirltrec(PIVector<float> *a, int astart, int astride, PIVector<float> *b, int bstart, int bstride, int m, int n);
+	void ftbase_ffttwcalc(PIVector<float> *a, int aoffset, int n1, int n2);
+};
+
+typedef PIFFT_double PIFFT;
+typedef PIFFT_double PIFFTd;
+typedef PIFFT_float PIFFTf;
+
 #endif // PIFFT_H