1889 lines
58 KiB
C++
1889 lines
58 KiB
C++
/*
|
|
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"
|
|
|
|
#ifndef MICRO_PIP
|
|
|
|
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);
|
|
assert(stackptr==0);
|
|
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);
|
|
assert(stackptr==0);
|
|
}
|
|
|
|
|
|
/*************************************************************************
|
|
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);
|
|
assertm(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan");
|
|
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);
|
|
assertm(stackptr==0, "Internal error in FTBaseGenerateComplexFFTPlan");
|
|
}
|
|
|
|
|
|
/*************************************************************************
|
|
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;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
#endif // MICRO_PIP
|