17template <
typename _Scalar>
20 typedef _Scalar Scalar;
21 typedef std::complex<Scalar> Complex;
22 std::vector<Complex> m_twiddles;
23 std::vector<int> m_stageRadix;
24 std::vector<int> m_stageRemainder;
25 std::vector<Complex> m_scratchBuf;
28 inline void make_twiddles(
int nfft,
bool inverse)
33 m_twiddles.resize(nfft);
34 double phinc = 0.25 * double(EIGEN_PI) / nfft;
35 Scalar flip =
inverse ? Scalar(1) : Scalar(-1);
36 m_twiddles[0] = Complex(Scalar(1), Scalar(0));
38 m_twiddles[nfft/2] = Complex(Scalar(-1), Scalar(0));
42 Scalar c = Scalar(
cos(i*8*phinc));
43 Scalar s = Scalar(
sin(i*8*phinc));
44 m_twiddles[i] = Complex(c, s*flip);
45 m_twiddles[nfft-i] = Complex(c, -s*flip);
49 Scalar c = Scalar(
cos((2*nfft-8*i)*phinc));
50 Scalar s = Scalar(
sin((2*nfft-8*i)*phinc));
51 m_twiddles[i] = Complex(s, c*flip);
52 m_twiddles[nfft-i] = Complex(s, -c*flip);
56 Scalar c = Scalar(
cos((8*i-2*nfft)*phinc));
57 Scalar s = Scalar(
sin((8*i-2*nfft)*phinc));
58 m_twiddles[i] = Complex(-s, c*flip);
59 m_twiddles[nfft-i] = Complex(-s, -c*flip);
63 Scalar c = Scalar(
cos((4*nfft-8*i)*phinc));
64 Scalar s = Scalar(
sin((4*nfft-8*i)*phinc));
65 m_twiddles[i] = Complex(-c, s*flip);
66 m_twiddles[nfft-i] = Complex(-c, -s*flip);
70 void factorize(
int nfft)
80 default: p += 2;
break;
86 m_stageRadix.push_back(p);
87 m_stageRemainder.push_back(n);
89 m_scratchBuf.resize(p);
93 template <
typename _Src>
95 void work(
int stage,Complex * xout,
const _Src * xin,
size_t fstride,
size_t in_stride)
97 int p = m_stageRadix[stage];
98 int m = m_stageRemainder[stage];
99 Complex * Fout_beg = xout;
100 Complex * Fout_end = xout + p*m;
108 work(stage+1, xout , xin, fstride*p,in_stride);
109 xin += fstride*in_stride;
110 }
while( (xout += m) != Fout_end );
114 xin += fstride*in_stride;
115 }
while(++xout != Fout_end );
121 case 2: bfly2(xout,fstride,m);
break;
122 case 3: bfly3(xout,fstride,m);
break;
123 case 4: bfly4(xout,fstride,m);
break;
124 case 5: bfly5(xout,fstride,m);
break;
125 default: bfly_generic(xout,fstride,m,p);
break;
130 void bfly2( Complex * Fout,
const size_t fstride,
int m)
132 for (
int k=0;k<m;++k) {
133 Complex t = Fout[m+k] * m_twiddles[k*fstride];
134 Fout[m+k] = Fout[k] - t;
140 void bfly4( Complex * Fout,
const size_t fstride,
const size_t m)
143 int negative_if_inverse = m_inverse * -2 +1;
144 for (
size_t k=0;k<m;++k) {
145 scratch[0] = Fout[k+m] * m_twiddles[k*fstride];
146 scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2];
147 scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3];
148 scratch[5] = Fout[k] - scratch[1];
150 Fout[k] += scratch[1];
151 scratch[3] = scratch[0] + scratch[2];
152 scratch[4] = scratch[0] - scratch[2];
153 scratch[4] = Complex( scratch[4].
imag()*negative_if_inverse , -scratch[4].
real()* negative_if_inverse );
155 Fout[k+2*m] = Fout[k] - scratch[3];
156 Fout[k] += scratch[3];
157 Fout[k+m] = scratch[5] + scratch[4];
158 Fout[k+3*m] = scratch[5] - scratch[4];
163 void bfly3( Complex * Fout,
const size_t fstride,
const size_t m)
166 const size_t m2 = 2*m;
170 epi3 = m_twiddles[fstride*m];
172 tw1=tw2=&m_twiddles[0];
175 scratch[1]=Fout[m] * *tw1;
176 scratch[2]=Fout[m2] * *tw2;
178 scratch[3]=scratch[1]+scratch[2];
179 scratch[0]=scratch[1]-scratch[2];
182 Fout[m] = Complex( Fout->real() - Scalar(.5)*scratch[3].
real() , Fout->imag() - Scalar(.5)*scratch[3].
imag() );
183 scratch[0] *= epi3.imag();
185 Fout[m2] = Complex( Fout[m].
real() + scratch[0].
imag() , Fout[m].
imag() - scratch[0].
real() );
186 Fout[m] += Complex( -scratch[0].
imag(),scratch[0].
real() );
192 void bfly5( Complex * Fout,
const size_t fstride,
const size_t m)
194 Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
197 Complex * twiddles = &m_twiddles[0];
200 ya = twiddles[fstride*m];
201 yb = twiddles[fstride*2*m];
210 for ( u=0; u<m; ++u ) {
213 scratch[1] = *Fout1 * tw[u*fstride];
214 scratch[2] = *Fout2 * tw[2*u*fstride];
215 scratch[3] = *Fout3 * tw[3*u*fstride];
216 scratch[4] = *Fout4 * tw[4*u*fstride];
218 scratch[7] = scratch[1] + scratch[4];
219 scratch[10] = scratch[1] - scratch[4];
220 scratch[8] = scratch[2] + scratch[3];
221 scratch[9] = scratch[2] - scratch[3];
223 *Fout0 += scratch[7];
224 *Fout0 += scratch[8];
226 scratch[5] = scratch[0] + Complex(
227 (scratch[7].
real()*ya.real() ) + (scratch[8].real() *yb.real() ),
228 (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real())
231 scratch[6] = Complex(
232 (scratch[10].
imag()*ya.imag()) + (scratch[9].imag()*yb.imag()),
233 -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag())
236 *Fout1 = scratch[5] - scratch[6];
237 *Fout4 = scratch[5] + scratch[6];
239 scratch[11] = scratch[0] +
241 (scratch[7].
real()*yb.real()) + (scratch[8].real()*ya.real()),
242 (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real())
245 scratch[12] = Complex(
246 -(scratch[10].
imag()*yb.imag()) + (scratch[9].imag()*ya.imag()),
247 (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag())
250 *Fout2=scratch[11]+scratch[12];
251 *Fout3=scratch[11]-scratch[12];
253 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
261 const size_t fstride,
267 Complex * twiddles = &m_twiddles[0];
269 int Norig =
static_cast<int>(m_twiddles.size());
270 Complex * scratchbuf = &m_scratchBuf[0];
272 for ( u=0; u<m; ++u ) {
274 for ( q1=0 ; q1<p ; ++q1 ) {
275 scratchbuf[q1] = Fout[ k ];
280 for ( q1=0 ; q1<p ; ++q1 ) {
282 Fout[ k ] = scratchbuf[0];
284 twidx +=
static_cast<int>(fstride) * k;
285 if (twidx>=Norig) twidx-=Norig;
286 t=scratchbuf[q] * twiddles[twidx];
295template <
typename _Scalar>
298 typedef _Scalar Scalar;
299 typedef std::complex<Scalar> Complex;
304 m_realTwiddles.clear();
308 void fwd( Complex * dst,
const Complex *src,
int nfft)
310 get_plan(nfft,
false).work(0, dst, src, 1,1);
314 void fwd2( Complex * dst,
const Complex *src,
int n0,
int n1)
316 EIGEN_UNUSED_VARIABLE(dst);
317 EIGEN_UNUSED_VARIABLE(src);
318 EIGEN_UNUSED_VARIABLE(n0);
319 EIGEN_UNUSED_VARIABLE(n1);
323 void inv2( Complex * dst,
const Complex *src,
int n0,
int n1)
325 EIGEN_UNUSED_VARIABLE(dst);
326 EIGEN_UNUSED_VARIABLE(src);
327 EIGEN_UNUSED_VARIABLE(n0);
328 EIGEN_UNUSED_VARIABLE(n1);
336 void fwd( Complex * dst,
const Scalar * src,
int nfft)
340 m_tmpBuf1.resize(nfft);
341 get_plan(nfft,
false).work(0, &m_tmpBuf1[0], src, 1,1);
342 std::copy(m_tmpBuf1.begin(),m_tmpBuf1.begin()+(nfft>>1)+1,dst );
345 int ncfft2 = nfft>>2;
346 Complex * rtw = real_twiddles(ncfft2);
349 fwd( dst,
reinterpret_cast<const Complex*
> (src), ncfft);
350 Complex dc(dst[0].
real() + dst[0].
imag());
351 Complex nyquist(dst[0].
real() - dst[0].
imag());
353 for ( k=1;k <= ncfft2 ; ++k ) {
354 Complex fpk = dst[k];
355 Complex fpnk =
conj(dst[ncfft-k]);
356 Complex f1k = fpk + fpnk;
357 Complex f2k = fpk - fpnk;
358 Complex tw= f2k * rtw[k-1];
359 dst[k] = (f1k + tw) * Scalar(.5);
360 dst[ncfft-k] =
conj(f1k -tw)*Scalar(.5);
363 dst[ncfft] = nyquist;
369 void inv(Complex * dst,
const Complex *src,
int nfft)
371 get_plan(nfft,
true).work(0, dst, src, 1,1);
376 void inv( Scalar * dst,
const Complex * src,
int nfft)
379 m_tmpBuf1.resize(nfft);
380 m_tmpBuf2.resize(nfft);
381 std::copy(src,src+(nfft>>1)+1,m_tmpBuf1.begin() );
382 for (
int k=1;k<(nfft>>1)+1;++k)
383 m_tmpBuf1[nfft-k] =
conj(m_tmpBuf1[k]);
384 inv(&m_tmpBuf2[0],&m_tmpBuf1[0],nfft);
385 for (
int k=0;k<nfft;++k)
386 dst[k] = m_tmpBuf2[k].
real();
390 int ncfft2 = nfft>>2;
391 Complex * rtw = real_twiddles(ncfft2);
392 m_tmpBuf1.resize(ncfft);
393 m_tmpBuf1[0] = Complex( src[0].
real() + src[ncfft].
real(), src[0].
real() - src[ncfft].
real() );
394 for (
int k = 1; k <= ncfft / 2; ++k) {
396 Complex fnkc =
conj(src[ncfft-k]);
397 Complex fek = fk + fnkc;
398 Complex tmp = fk - fnkc;
399 Complex fok = tmp *
conj(rtw[k-1]);
400 m_tmpBuf1[k] = fek + fok;
401 m_tmpBuf1[ncfft-k] =
conj(fek - fok);
403 get_plan(ncfft,
true).work(0,
reinterpret_cast<Complex*
>(dst), &m_tmpBuf1[0], 1,1);
408 typedef kiss_cpx_fft<Scalar> PlanData;
409 typedef std::map<int,PlanData> PlanMap;
412 std::map<int, std::vector<Complex> > m_realTwiddles;
413 std::vector<Complex> m_tmpBuf1;
414 std::vector<Complex> m_tmpBuf2;
417 int PlanKey(
int nfft,
bool isinverse)
const {
return (nfft<<1) | int(isinverse); }
420 PlanData & get_plan(
int nfft,
bool inverse)
423 PlanData & pd = m_plans[ PlanKey(nfft,
inverse) ];
424 if ( pd.m_twiddles.size() == 0 ) {
425 pd.make_twiddles(nfft,
inverse);
432 Complex * real_twiddles(
int ncfft2)
435 std::vector<Complex> & twidref = m_realTwiddles[ncfft2];
436 if ( (
int)twidref.size() != ncfft2 ) {
437 twidref.resize(ncfft2);
438 int ncfft= ncfft2<<1;
439 Scalar pi =
acos( Scalar(-1) );
440 for (
int k=1;k<=ncfft2;++k)
441 twidref[k-1] =
exp( Complex(0,-pi * (Scalar(k) / ncfft + Scalar(.5)) ) );
Namespace containing all symbols from the Eigen library.
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_cos_op< typename Derived::Scalar >, const Derived > cos(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_inverse_op< typename Derived::Scalar >, const Derived > inverse(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_exp_op< typename Derived::Scalar >, const Derived > exp(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_imag_op< typename Derived::Scalar >, const Derived > imag(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sin_op< typename Derived::Scalar >, const Derived > sin(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_acos_op< typename Derived::Scalar >, const Derived > acos(const Eigen::ArrayBase< Derived > &x)