tree changes
This commit is contained in:
1893
libs/main/math/pifft.cpp
Normal file
1893
libs/main/math/pifft.cpp
Normal file
@@ -0,0 +1,1893 @@
|
||||
/*
|
||||
PIP - Platform Independent Primitives
|
||||
Class for FFT, IFFT and Hilbert transformations
|
||||
Andrey Bychkov work.a.b@yandex.ru, Ivan Pelipenko peri4ko@yandex.ru
|
||||
|
||||
This program is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU Lesser General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "pifft.h"
|
||||
|
||||
|
||||
PIFFT_double::PIFFT_double() {
|
||||
}
|
||||
|
||||
|
||||
PIVector<complexd> * PIFFT_double::calcFFT(const PIVector<complexd> & val) {
|
||||
result.clear();
|
||||
if (val.size_s() < 4) return &result;
|
||||
fftc1d(val, val.size());
|
||||
return &result;
|
||||
}
|
||||
|
||||
|
||||
PIVector<complexd> * PIFFT_double::calcFFTinverse(const PIVector<complexd> & val) {
|
||||
result.clear();
|
||||
if (val.size_s() < 4) return &result;
|
||||
fftc1dinv(val, val.size());
|
||||
return &result;
|
||||
}
|
||||
|
||||
|
||||
PIVector<complexd> * PIFFT_double::calcHilbert(const PIVector<double> & 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.;
|
||||
for (uint i = result.size() / 2; i < result.size(); i++) result[i] = 0;
|
||||
fftc1dinv(result, result.size());
|
||||
return &result;
|
||||
}
|
||||
|
||||
|
||||
PIVector<complexd> * PIFFT_double::calcFFT(const PIVector<double> & val) {
|
||||
result.clear();
|
||||
if (val.size_s() < 4) return &result;
|
||||
fftc1r(val, val.size());
|
||||
return &result;
|
||||
}
|
||||
|
||||
|
||||
PIVector<double> PIFFT_double::getAmplitude() const {
|
||||
PIVector<double> ret;
|
||||
ret.resize(result.size());
|
||||
double tmp;
|
||||
for (uint i = 0; i < result.size(); i++) {
|
||||
tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
|
||||
ret[i] = tmp;
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
PIVector<double> PIFFT_double::getReal() const {
|
||||
PIVector<double> ret;
|
||||
ret.resize(result.size());
|
||||
for (uint i = 0; i < result.size(); i++)
|
||||
ret[i] = result[i].real();
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
PIVector<double> PIFFT_double::getImag() const {
|
||||
PIVector<double> ret;
|
||||
ret.resize(result.size());
|
||||
for (uint i = 0; i < result.size(); i++)
|
||||
ret[i] = result[i].imag();
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
void PIFFT_double::fftc1d(const PIVector<complexd> & a, uint n) {
|
||||
createPlan(n);
|
||||
uint i;
|
||||
PIVector<double> 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] = complexd(buf[2 * i + 0], buf[2 * i + 1]);
|
||||
}
|
||||
|
||||
|
||||
void PIFFT_double::fftc1r(const PIVector<double> & a, uint n) {
|
||||
uint i;
|
||||
if (n % 2 == 0) {
|
||||
PIVector<double> buf;
|
||||
uint n2 = n / 2;
|
||||
buf = a;
|
||||
createPlan(n2);
|
||||
ftbaseexecuteplan(&buf, 0, n2, &curplan);
|
||||
result.resize(n);
|
||||
uint idx;
|
||||
complexd hn, hmnc, v;
|
||||
for (i = 0; i <= n2; i++) {
|
||||
idx = 2 * (i % n2);
|
||||
hn = complexd(buf[idx + 0], buf[idx + 1]);
|
||||
idx = 2 * ((n2 - i) % n2);
|
||||
hmnc = complexd(buf[idx + 0], -buf[idx + 1]);
|
||||
v = complexd(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<complexd> cbuf;
|
||||
cbuf.resize(n);
|
||||
for (i = 0; i < n; i++)
|
||||
cbuf[i] = complexd(a[i], 0.);
|
||||
fftc1d(cbuf, n);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void PIFFT_double::fftc1dinv(const PIVector<complexd> & a, uint n) {
|
||||
PIVector<complexd> 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] / (double)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);
|
||||
}
|
||||
|
||||
|
||||
void PIFFT_double::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_double::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 = 8;
|
||||
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(8 * (*planarraysize));
|
||||
*planarraysize = 8 * (*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_double::ftbase_ftbaseprecomputeplanrec(ftplan * plan,
|
||||
int entryoffset,
|
||||
ae_int_t stackptr) {
|
||||
int n1, n2, n, m, offs;
|
||||
double 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_double::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_double::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_double::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_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_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)
|
||||
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_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_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;
|
||||
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_double::ftbaseexecuteplan(PIVector<double> * 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 doubles)
|
||||
|
||||
-- ALGLIB --
|
||||
Copyright 01.05.2009 by Bochkanov Sergey
|
||||
*************************************************************************/
|
||||
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;
|
||||
double t1x, t1y, t2x, t2y, t3x, t3y, t4x, t4y, t5x, t5y;
|
||||
double m1x, m1y, m2x, m2y, m3x, m3y, m4x, m4y, m5x, m5y;
|
||||
double s1x, s1y, s2x, s2y, s3x, s3y, s4x, s4y, s5x, s5y;
|
||||
double 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<double> & 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_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;
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//=================================================================================================
|
||||
//*************************************************************************************************
|
||||
//---------------------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() const {
|
||||
PIVector<float> ret;
|
||||
ret.resize(result.size());
|
||||
float tmp;
|
||||
for (uint i = 0; i < result.size(); i++) {
|
||||
tmp = sqrt(result[i].real() * result[i].real() + result[i].imag() * result[i].imag());
|
||||
ret[i] = tmp;
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
PIVector<float> PIFFT_float::getReal() const {
|
||||
PIVector<float> ret;
|
||||
ret.resize(result.size());
|
||||
for (uint i = 0; i < result.size(); i++)
|
||||
ret[i] = result[i].real();
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
PIVector<float> PIFFT_float::getImag() const {
|
||||
PIVector<float> ret;
|
||||
ret.resize(result.size());
|
||||
for (uint i = 0; i < result.size(); i++)
|
||||
ret[i] = result[i].imag();
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user