From f0d4d75f651fa547e0554e1d8fe6c25fae025097 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=91=D1=8B=D1=87=D0=BA=D0=BE=D0=B2=20=D0=90=D0=BD=D0=B4?= =?UTF-8?q?=D1=80=D0=B5=D0=B9?= Date: Mon, 4 May 2015 09:54:02 +0000 Subject: [PATCH] PIFFT_float git-svn-id: svn://db.shs.com.ru/pip@124 12ceb7fc-bf1f-11e4-8940-5bc7170c53b5 --- src/math/pifft.cpp | 970 +++++++++++++++++++++++++++++++++++++++++++-- src/math/pifft.h | 57 ++- 2 files changed, 995 insertions(+), 32 deletions(-) diff --git a/src/math/pifft.cpp b/src/math/pifft.cpp index cf6df93a..ef5807bd 100644 --- 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 * PIFFT::calcFFT(const PIVector & val) { +PIVector * PIFFT_double::calcFFT(const PIVector & val) { result.clear(); if (val.size_s() < 4) return &result; fftc1d(val, val.size()); @@ -33,7 +32,7 @@ PIVector * PIFFT::calcFFT(const PIVector & val) { } -PIVector * PIFFT::calcFFTinverse(const PIVector & val) { +PIVector * PIFFT_double::calcFFTinverse(const PIVector & val) { result.clear(); if (val.size_s() < 4) return &result; fftc1dinv(val, val.size()); @@ -41,7 +40,7 @@ PIVector * PIFFT::calcFFTinverse(const PIVector & val) { } -PIVector * PIFFT::calcHilbert(const PIVector & val) { +PIVector * PIFFT_double::calcHilbert(const PIVector & val) { result.clear(); if (val.size_s() < 4) return &result; fftc1r(val, val.size()); @@ -52,7 +51,7 @@ PIVector * PIFFT::calcHilbert(const PIVector & val) { } -PIVector * PIFFT::calcFFT(const PIVector & val) { +PIVector * PIFFT_double::calcFFT(const PIVector & val) { result.clear(); if (val.size_s() < 4) return &result; fftc1r(val, val.size()); @@ -60,7 +59,7 @@ PIVector * PIFFT::calcFFT(const PIVector & val) { } -PIVector PIFFT::getAmplitude() { +PIVector PIFFT_double::getAmplitude() { PIVector a; double tmp; for (uint i = 0; i < result.size(); i++) { @@ -71,7 +70,7 @@ PIVector PIFFT::getAmplitude() { } -void PIFFT::fftc1d(const PIVector & a, uint n) { +void PIFFT_double::fftc1d(const PIVector & a, uint n) { createPlan(n); uint i; PIVector buf; @@ -87,7 +86,7 @@ void PIFFT::fftc1d(const PIVector & a, uint n) { } -void PIFFT::fftc1r(const PIVector & a, uint n) { +void PIFFT_double::fftc1r(const PIVector & a, uint n) { uint i; if (n % 2 == 0) { PIVector buf; @@ -119,7 +118,7 @@ void PIFFT::fftc1r(const PIVector & a, uint n) { } -void PIFFT::fftc1dinv(const PIVector & a, uint n) { +void PIFFT_double::fftc1dinv(const PIVector & a, uint n) { PIVector cbuf; cbuf.resize(n); uint i; @@ -133,18 +132,17 @@ void PIFFT::fftc1dinv(const PIVector & 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 * a, int m, int n, int astart, PIVector * buf) { +void PIFFT_double::ftbase_internalreallintranspose(PIVector * a, int m, int n, int astart, PIVector * 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 * a, int astart, int astride, PIVector * b, int bstart, int bstride, int m, int n) { +void PIFFT_double::ftbase_fftirltrec(PIVector * a, int astart, int astride, PIVector * 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 * a, int astart, int astride, PIV } -void PIFFT::ftbase_internalcomplexlintranspose(PIVector * a, int m, int n, int astart, PIVector * buf) { +void PIFFT_double::ftbase_internalcomplexlintranspose(PIVector * a, int m, int n, int astart, PIVector * 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 * a, int astart, int astride, PIVector * b, int bstart, int bstride, int m, int n) { +void PIFFT_double::ftbase_ffticltrec(PIVector * a, int astart, int astride, PIVector * 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 * a, int astart, int astride, PIV } -void PIFFT::ftbaseexecuteplan(PIVector * a, int aoffset, int n, ftplan * plan) { +void PIFFT_double::ftbaseexecuteplan(PIVector * 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 * a, int aoffset, ftplan * plan, int entryoffset, ae_int_t stackptr) { +void PIFFT_double::ftbaseexecuteplanrec(PIVector * 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 * a, int aoffset, int n1, int n2) { +void PIFFT_double::ftbase_ffttwcalc(PIVector * 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 * a, int aoffset, int n1, int n2) } } } + + +//================================================================================================= +//************************************************************************************************* +//---------------------FLOAT-----------------------FLOAT------------------------------------------- +//************************************************************************************************* +//================================================================================================= + + +PIFFT_float::PIFFT_float() { +} + + +PIVector * PIFFT_float::calcFFT(const PIVector & val) { + result.clear(); + if (val.size_s() < 4) return &result; + fftc1d(val, val.size()); + return &result; +} + + +PIVector * PIFFT_float::calcFFTinverse(const PIVector & val) { + result.clear(); + if (val.size_s() < 4) return &result; + fftc1dinv(val, val.size()); + return &result; +} + + +PIVector * PIFFT_float::calcHilbert(const PIVector & 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 * PIFFT_float::calcFFT(const PIVector & val) { + result.clear(); + if (val.size_s() < 4) return &result; + fftc1r(val, val.size()); + return &result; +} + + +PIVector PIFFT_float::getAmplitude() { + PIVector 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 & a, uint n) { + createPlan(n); + uint i; + PIVector 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 & a, uint n) { + uint i; + if (n % 2 == 0) { + PIVector 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 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 & a, uint n) { + PIVector 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(stackmemsize, 1)); //ae_vector_set_length(&curplan.stackbuf, ae_maxint(stackmemsize, 1)); + curplan.tmpbuf.resize(piMax(tmpmemsize, 1)); //ae_vector_set_length(&(curplan.tmpbuf), ae_maxint(tmpmemsize, 1)); + curplan.precomputed.resize(piMax(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(*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(*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(*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 * a, int m, int n, int astart, PIVector * 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 * a, int astart, int astride, PIVector * 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 * a, int m, int n, int astart, PIVector * 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 * a, int astart, int astride, PIVector * b, int bstart, int bstride, int m, int n) { + int idx1, idx2, m2, m1, n1; + if (m == 0 || n == 0) + return; + if (piMax(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 * 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 * 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 & 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 * 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; + } + } + } +} + diff --git a/src/math/pifft.h b/src/math/pifft.h index bb45b4bf..73c75a41 100644 --- 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 * calcFFT(const PIVector &val); PIVector * calcFFT(const PIVector &val); @@ -38,11 +38,7 @@ public: private: PIVector 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 plan; @@ -71,7 +67,54 @@ private: void ftbase_internalreallintranspose(PIVector *a, int m, int n, int astart, PIVector *buf); void ftbase_fftirltrec(PIVector *a, int astart, int astride, PIVector *b, int bstart, int bstride, int m, int n); void ftbase_ffttwcalc(PIVector *a, int aoffset, int n1, int n2); - }; +class PIP_EXPORT PIFFT_float +{ +public: + PIFFT_float(); + + PIVector * calcFFT(const PIVector &val); + PIVector * calcFFT(const PIVector &val); + PIVector * calcFFTinverse(const PIVector &val); + PIVector * calcHilbert(const PIVector &val); + PIVector getAmplitude(); + +private: + PIVector result; + typedef ptrdiff_t ae_int_t; + + struct ftplan { + PIVector plan; + PIVector precomputed; + PIVector tmpbuf; + PIVector stackbuf; + }; + + ftplan curplan; + + void fftc1d(const PIVector &a, uint n); + void fftc1r(const PIVector &a, uint n); + void fftc1dinv(const PIVector &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 *a, int aoffset, int n, ftplan *plan); + void ftbaseexecuteplanrec(PIVector *a, int aoffset, ftplan *plan, int entryoffset, ae_int_t stackptr); + void ftbase_internalcomplexlintranspose(PIVector *a, int m, int n, int astart, PIVector *buf); + void ftbase_ffticltrec(PIVector *a, int astart, int astride, PIVector *b, int bstart, int bstride, int m, int n); + void ftbase_internalreallintranspose(PIVector *a, int m, int n, int astart, PIVector *buf); + void ftbase_fftirltrec(PIVector *a, int astart, int astride, PIVector *b, int bstart, int bstride, int m, int n); + void ftbase_ffttwcalc(PIVector *a, int aoffset, int n1, int n2); +}; + +typedef PIFFT_double PIFFT; +typedef PIFFT_double PIFFTd; +typedef PIFFT_float PIFFTf; + #endif // PIFFT_H