Commit ef0739a2 authored by Martin Reinecke's avatar Martin Reinecke

reduce memory consumption of Bluestein plans

parent c0f74f61
......@@ -2498,8 +2498,14 @@ template<typename T0> class fftblue
plan.exec (akf.data(),1.,true);
/* do the convolution */
for (size_t m=0; m<n2; ++m)
akf[0] = akf[0].template special_mul<!fwd>(bkf[0]);
for (size_t m=1; m<(n2+1)/2; ++m)
{
akf[m] = akf[m].template special_mul<!fwd>(bkf[m]);
akf[n2-m] = akf[n2-m].template special_mul<!fwd>(bkf[m]);
}
if ((n2&1)==0)
akf[n2/2] = akf[n2/2].template special_mul<!fwd>(bkf[n2/2]);
/* inverse FFT */
plan.exec (akf.data(),1.,false);
......@@ -2511,7 +2517,7 @@ template<typename T0> class fftblue
public:
POCKETFFT_NOINLINE fftblue(size_t length)
: n(length), n2(util::good_size_cmplx(n*2-1)), plan(n2), mem(n+n2),
: n(length), n2(util::good_size_cmplx(n*2-1)), plan(n2), mem(n+n2/2+1),
bk(mem.data()), bkf(mem.data()+n)
{
/* initialize b_k */
......@@ -2528,13 +2534,16 @@ template<typename T0> class fftblue
}
/* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */
arr<cmplx<T0>> tbkf(n2);
T0 xn2 = T0(1)/T0(n2);
bkf[0] = bk[0]*xn2;
tbkf[0] = bk[0]*xn2;
for (size_t m=1; m<n; ++m)
bkf[m] = bkf[n2-m] = bk[m]*xn2;
tbkf[m] = tbkf[n2-m] = bk[m]*xn2;
for (size_t m=n;m<=(n2-n);++m)
bkf[m].Set(0.,0.);
plan.exec(bkf,1.,true);
tbkf[m].Set(0.,0.);
plan.exec(tbkf.data(),1.,true);
for (size_t i=0; i<n2/2+1; ++i)
bkf[i] = tbkf[i];
}
template<typename T> void exec(cmplx<T> c[], T0 fct, bool fwd) const
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment