diff options
author | Mark Borgerding <mark@borgerding.net> | 2009-10-31 00:13:22 -0400 |
---|---|---|
committer | Mark Borgerding <mark@borgerding.net> | 2009-10-31 00:13:22 -0400 |
commit | ec70f8006be1e46055a0d65850c8c60c3bebf6ed (patch) | |
tree | 32621c8a9a67827165cd02e3caefa45347c8f411 /unsupported | |
parent | 4c3345364e079429dcfc17da63364ee75b9c0636 (diff) |
added inlines to a bunch of functions
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/Eigen/FFT | 12 | ||||
-rw-r--r-- | unsupported/Eigen/src/FFT/ei_fftw_impl.h | 3 | ||||
-rw-r--r-- | unsupported/Eigen/src/FFT/ei_kissfft_impl.h | 695 |
3 files changed, 365 insertions, 345 deletions
diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index b1d3d9f0e..8f7a2fcae 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -85,6 +85,7 @@ class FFT inline void ClearFlag(Flag f) { m_flag &= (~(int)f);} + inline void fwd( Complex * dst, const Scalar * src, int nfft) { m_impl.fwd(dst,src,nfft); @@ -92,12 +93,14 @@ class FFT ReflectSpectrum(dst,nfft); } + inline void fwd( Complex * dst, const Complex * src, int nfft) { m_impl.fwd(dst,src,nfft); } template <typename _Input> + inline void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src) { if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) ) @@ -108,6 +111,7 @@ class FFT } template<typename InputDerived, typename ComplexDerived> + inline void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived) @@ -125,6 +129,7 @@ class FFT fwd( &dst[0],&src[0],src.size() ); } + inline void inv( Complex * dst, const Complex * src, int nfft) { m_impl.inv( dst,src,nfft ); @@ -132,6 +137,7 @@ class FFT scale(dst,1./nfft,nfft); } + inline void inv( Scalar * dst, const Complex * src, int nfft) { m_impl.inv( dst,src,nfft ); @@ -140,6 +146,7 @@ class FFT } template<typename OutputDerived, typename ComplexDerived> + inline void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived) @@ -157,6 +164,7 @@ class FFT } template <typename _Output> + inline void inv( std::vector<_Output> & dst, const std::vector<Complex> & src) { if ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) @@ -171,18 +179,22 @@ class FFT // TODO: handle Eigen MatrixBase // ---> i added fwd and inv specializations above + unit test, is this enough? (bjacob) + inline impl_type & impl() {return m_impl;} private: template <typename _It,typename _Val> + inline void scale(_It x,_Val s,int nx) { for (int k=0;k<nx;++k) *x++ *= s; } + inline void ReflectSpectrum(Complex * freq,int nfft) { + // create the implicit right-half spectrum (conjugate-mirror of the left-half) int nhbins=(nfft>>1)+1; for (int k=nhbins;k < nfft; ++k ) freq[k] = conj(freq[nfft-k]); diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h index 18473a29b..a66b7398c 100644 --- a/unsupported/Eigen/src/FFT/ei_fftw_impl.h +++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h @@ -166,6 +166,7 @@ m_plans.clear(); } + // complex-to-complex forward FFT inline void fwd( Complex * dst,const Complex *src,int nfft) { @@ -208,3 +209,5 @@ return m_plans[key]; } }; +/* vim: set filetype=cpp et sw=2 ts=2 ai: */ + diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h index 091e730d1..5c958d1ec 100644 --- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h +++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h @@ -27,379 +27,384 @@ // This FFT implementation was derived from kissfft http:sourceforge.net/projects/kissfft // Copyright 2003-2009 Mark Borgerding - template <typename _Scalar> - struct ei_kiss_cpx_fft +template <typename _Scalar> +struct ei_kiss_cpx_fft +{ + typedef _Scalar Scalar; + typedef std::complex<Scalar> Complex; + std::vector<Complex> m_twiddles; + std::vector<int> m_stageRadix; + std::vector<int> m_stageRemainder; + std::vector<Complex> m_scratchBuf; + bool m_inverse; + + inline + void make_twiddles(int nfft,bool inverse) { - typedef _Scalar Scalar; - typedef std::complex<Scalar> Complex; - std::vector<Complex> m_twiddles; - std::vector<int> m_stageRadix; - std::vector<int> m_stageRemainder; - std::vector<Complex> m_scratchBuf; - bool m_inverse; - - void make_twiddles(int nfft,bool inverse) - { - m_inverse = inverse; - m_twiddles.resize(nfft); - Scalar phinc = (inverse?2:-2)* acos( (Scalar) -1) / nfft; - for (int i=0;i<nfft;++i) - m_twiddles[i] = exp( Complex(0,i*phinc) ); + m_inverse = inverse; + m_twiddles.resize(nfft); + Scalar phinc = (inverse?2:-2)* acos( (Scalar) -1) / nfft; + for (int i=0;i<nfft;++i) + m_twiddles[i] = exp( Complex(0,i*phinc) ); + } + + void factorize(int nfft) + { + //start factoring out 4's, then 2's, then 3,5,7,9,... + int n= nfft; + int p=4; + do { + while (n % p) { + switch (p) { + case 4: p = 2; break; + case 2: p = 3; break; + default: p += 2; break; + } + if (p*p>n) + p=n;// impossible to have a factor > sqrt(n) } + n /= p; + m_stageRadix.push_back(p); + m_stageRemainder.push_back(n); + if ( p > 5 ) + m_scratchBuf.resize(p); // scratchbuf will be needed in bfly_generic + }while(n>1); + } + + template <typename _Src> + inline + void work( int stage,Complex * xout, const _Src * xin, size_t fstride,size_t in_stride) + { + int p = m_stageRadix[stage]; + int m = m_stageRemainder[stage]; + Complex * Fout_beg = xout; + Complex * Fout_end = xout + p*m; - void factorize(int nfft) - { - //start factoring out 4's, then 2's, then 3,5,7,9,... - int n= nfft; - int p=4; - do { - while (n % p) { - switch (p) { - case 4: p = 2; break; - case 2: p = 3; break; - default: p += 2; break; - } - if (p*p>n) - p=n;// impossible to have a factor > sqrt(n) - } - n /= p; - m_stageRadix.push_back(p); - m_stageRemainder.push_back(n); - if ( p > 5 ) - m_scratchBuf.resize(p); // scratchbuf will be needed in bfly_generic - }while(n>1); + if (m>1) { + do{ + // recursive call: + // DFT of size m*p performed by doing + // p instances of smaller DFTs of size m, + // each one takes a decimated version of the input + work(stage+1, xout , xin, fstride*p,in_stride); + xin += fstride*in_stride; + }while( (xout += m) != Fout_end ); + }else{ + do{ + *xout = *xin; + xin += fstride*in_stride; + }while(++xout != Fout_end ); } - - template <typename _Src> - void work( int stage,Complex * xout, const _Src * xin, size_t fstride,size_t in_stride) - { - int p = m_stageRadix[stage]; - int m = m_stageRemainder[stage]; - Complex * Fout_beg = xout; - Complex * Fout_end = xout + p*m; - - if (m>1) { - do{ - // recursive call: - // DFT of size m*p performed by doing - // p instances of smaller DFTs of size m, - // each one takes a decimated version of the input - work(stage+1, xout , xin, fstride*p,in_stride); - xin += fstride*in_stride; - }while( (xout += m) != Fout_end ); - }else{ - do{ - *xout = *xin; - xin += fstride*in_stride; - }while(++xout != Fout_end ); - } - xout=Fout_beg; - - // recombine the p smaller DFTs - switch (p) { - case 2: bfly2(xout,fstride,m); break; - case 3: bfly3(xout,fstride,m); break; - case 4: bfly4(xout,fstride,m); break; - case 5: bfly5(xout,fstride,m); break; - default: bfly_generic(xout,fstride,m,p); break; - } - } - - inline - void bfly2( Complex * Fout, const size_t fstride, int m) - { - for (int k=0;k<m;++k) { - Complex t = Fout[m+k] * m_twiddles[k*fstride]; - Fout[m+k] = Fout[k] - t; - Fout[k] += t; - } + xout=Fout_beg; + + // recombine the p smaller DFTs + switch (p) { + case 2: bfly2(xout,fstride,m); break; + case 3: bfly3(xout,fstride,m); break; + case 4: bfly4(xout,fstride,m); break; + case 5: bfly5(xout,fstride,m); break; + default: bfly_generic(xout,fstride,m,p); break; } + } - inline - void bfly4( Complex * Fout, const size_t fstride, const size_t m) - { - Complex scratch[6]; - int negative_if_inverse = m_inverse * -2 +1; - for (size_t k=0;k<m;++k) { - scratch[0] = Fout[k+m] * m_twiddles[k*fstride]; - scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2]; - scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3]; - scratch[5] = Fout[k] - scratch[1]; - - Fout[k] += scratch[1]; - scratch[3] = scratch[0] + scratch[2]; - scratch[4] = scratch[0] - scratch[2]; - scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse ); - - Fout[k+2*m] = Fout[k] - scratch[3]; - Fout[k] += scratch[3]; - Fout[k+m] = scratch[5] + scratch[4]; - Fout[k+3*m] = scratch[5] - scratch[4]; - } + inline + void bfly2( Complex * Fout, const size_t fstride, int m) + { + for (int k=0;k<m;++k) { + Complex t = Fout[m+k] * m_twiddles[k*fstride]; + Fout[m+k] = Fout[k] - t; + Fout[k] += t; } + } - inline - void bfly3( Complex * Fout, const size_t fstride, const size_t m) - { - size_t k=m; - const size_t m2 = 2*m; - Complex *tw1,*tw2; - Complex scratch[5]; - Complex epi3; - epi3 = m_twiddles[fstride*m]; - - tw1=tw2=&m_twiddles[0]; - - do{ - scratch[1]=Fout[m] * *tw1; - scratch[2]=Fout[m2] * *tw2; - - scratch[3]=scratch[1]+scratch[2]; - scratch[0]=scratch[1]-scratch[2]; - tw1 += fstride; - tw2 += fstride*2; - Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() ); - scratch[0] *= epi3.imag(); - *Fout += scratch[3]; - Fout[m2] = Complex( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() ); - Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() ); - ++Fout; - }while(--k); + inline + void bfly4( Complex * Fout, const size_t fstride, const size_t m) + { + Complex scratch[6]; + int negative_if_inverse = m_inverse * -2 +1; + for (size_t k=0;k<m;++k) { + scratch[0] = Fout[k+m] * m_twiddles[k*fstride]; + scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2]; + scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3]; + scratch[5] = Fout[k] - scratch[1]; + + Fout[k] += scratch[1]; + scratch[3] = scratch[0] + scratch[2]; + scratch[4] = scratch[0] - scratch[2]; + scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse ); + + Fout[k+2*m] = Fout[k] - scratch[3]; + Fout[k] += scratch[3]; + Fout[k+m] = scratch[5] + scratch[4]; + Fout[k+3*m] = scratch[5] - scratch[4]; } + } - inline - void bfly5( Complex * Fout, const size_t fstride, const size_t m) - { - Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; - size_t u; - Complex scratch[13]; - Complex * twiddles = &m_twiddles[0]; - Complex *tw; - Complex ya,yb; - ya = twiddles[fstride*m]; - yb = twiddles[fstride*2*m]; - - Fout0=Fout; - Fout1=Fout0+m; - Fout2=Fout0+2*m; - Fout3=Fout0+3*m; - Fout4=Fout0+4*m; - - tw=twiddles; - for ( u=0; u<m; ++u ) { - scratch[0] = *Fout0; - - scratch[1] = *Fout1 * tw[u*fstride]; - scratch[2] = *Fout2 * tw[2*u*fstride]; - scratch[3] = *Fout3 * tw[3*u*fstride]; - scratch[4] = *Fout4 * tw[4*u*fstride]; - - scratch[7] = scratch[1] + scratch[4]; - scratch[10] = scratch[1] - scratch[4]; - scratch[8] = scratch[2] + scratch[3]; - scratch[9] = scratch[2] - scratch[3]; - - *Fout0 += scratch[7]; - *Fout0 += scratch[8]; - - scratch[5] = scratch[0] + Complex( - (scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ), - (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real()) - ); - - scratch[6] = Complex( - (scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()), - -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag()) + inline + void bfly3( Complex * Fout, const size_t fstride, const size_t m) + { + size_t k=m; + const size_t m2 = 2*m; + Complex *tw1,*tw2; + Complex scratch[5]; + Complex epi3; + epi3 = m_twiddles[fstride*m]; + + tw1=tw2=&m_twiddles[0]; + + do{ + scratch[1]=Fout[m] * *tw1; + scratch[2]=Fout[m2] * *tw2; + + scratch[3]=scratch[1]+scratch[2]; + scratch[0]=scratch[1]-scratch[2]; + tw1 += fstride; + tw2 += fstride*2; + Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() ); + scratch[0] *= epi3.imag(); + *Fout += scratch[3]; + Fout[m2] = Complex( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() ); + Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() ); + ++Fout; + }while(--k); + } + + inline + void bfly5( Complex * Fout, const size_t fstride, const size_t m) + { + Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; + size_t u; + Complex scratch[13]; + Complex * twiddles = &m_twiddles[0]; + Complex *tw; + Complex ya,yb; + ya = twiddles[fstride*m]; + yb = twiddles[fstride*2*m]; + + Fout0=Fout; + Fout1=Fout0+m; + Fout2=Fout0+2*m; + Fout3=Fout0+3*m; + Fout4=Fout0+4*m; + + tw=twiddles; + for ( u=0; u<m; ++u ) { + scratch[0] = *Fout0; + + scratch[1] = *Fout1 * tw[u*fstride]; + scratch[2] = *Fout2 * tw[2*u*fstride]; + scratch[3] = *Fout3 * tw[3*u*fstride]; + scratch[4] = *Fout4 * tw[4*u*fstride]; + + scratch[7] = scratch[1] + scratch[4]; + scratch[10] = scratch[1] - scratch[4]; + scratch[8] = scratch[2] + scratch[3]; + scratch[9] = scratch[2] - scratch[3]; + + *Fout0 += scratch[7]; + *Fout0 += scratch[8]; + + scratch[5] = scratch[0] + Complex( + (scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ), + (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real()) + ); + + scratch[6] = Complex( + (scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()), + -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag()) + ); + + *Fout1 = scratch[5] - scratch[6]; + *Fout4 = scratch[5] + scratch[6]; + + scratch[11] = scratch[0] + + Complex( + (scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()), + (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real()) ); - *Fout1 = scratch[5] - scratch[6]; - *Fout4 = scratch[5] + scratch[6]; - - scratch[11] = scratch[0] + - Complex( - (scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()), - (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real()) - ); - - scratch[12] = Complex( - -(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()), - (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag()) - ); + scratch[12] = Complex( + -(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()), + (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag()) + ); - *Fout2=scratch[11]+scratch[12]; - *Fout3=scratch[11]-scratch[12]; + *Fout2=scratch[11]+scratch[12]; + *Fout3=scratch[11]-scratch[12]; - ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; - } + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; } + } + + /* perform the butterfly for one stage of a mixed radix FFT */ + inline + void bfly_generic( + Complex * Fout, + const size_t fstride, + int m, + int p + ) + { + int u,k,q1,q; + Complex * twiddles = &m_twiddles[0]; + Complex t; + int Norig = m_twiddles.size(); + Complex * scratchbuf = &m_scratchBuf[0]; + + for ( u=0; u<m; ++u ) { + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + scratchbuf[q1] = Fout[ k ]; + k += m; + } - /* perform the butterfly for one stage of a mixed radix FFT */ - inline - void bfly_generic( - Complex * Fout, - const size_t fstride, - int m, - int p - ) - { - int u,k,q1,q; - Complex * twiddles = &m_twiddles[0]; - Complex t; - int Norig = m_twiddles.size(); - Complex * scratchbuf = &m_scratchBuf[0]; - - for ( u=0; u<m; ++u ) { - k=u; - for ( q1=0 ; q1<p ; ++q1 ) { - scratchbuf[q1] = Fout[ k ]; - k += m; - } - - k=u; - for ( q1=0 ; q1<p ; ++q1 ) { - int twidx=0; - Fout[ k ] = scratchbuf[0]; - for (q=1;q<p;++q ) { - twidx += fstride * k; - if (twidx>=Norig) twidx-=Norig; - t=scratchbuf[q] * twiddles[twidx]; - Fout[ k ] += t; - } - k += m; + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + int twidx=0; + Fout[ k ] = scratchbuf[0]; + for (q=1;q<p;++q ) { + twidx += fstride * k; + if (twidx>=Norig) twidx-=Norig; + t=scratchbuf[q] * twiddles[twidx]; + Fout[ k ] += t; } + k += m; } } - }; - - template <typename _Scalar> - struct ei_kissfft_impl + } +}; + +template <typename _Scalar> +struct ei_kissfft_impl +{ + typedef _Scalar Scalar; + typedef std::complex<Scalar> Complex; + + void clear() + { + m_plans.clear(); + m_realTwiddles.clear(); + } + + inline + void fwd( Complex * dst,const Complex *src,int nfft) { - typedef _Scalar Scalar; - typedef std::complex<Scalar> Complex; - - void clear() - { - m_plans.clear(); - m_realTwiddles.clear(); + get_plan(nfft,false).work(0, dst, src, 1,1); + } + + // real-to-complex forward FFT + // perform two FFTs of src even and src odd + // then twiddle to recombine them into the half-spectrum format + // then fill in the conjugate symmetric half + inline + void fwd( Complex * dst,const Scalar * src,int nfft) + { + if ( nfft&3 ) { + // use generic mode for odd + m_tmpBuf1.resize(nfft); + get_plan(nfft,false).work(0, &m_tmpBuf1[0], src, 1,1); + std::copy(m_tmpBuf1.begin(),m_tmpBuf1.begin()+(nfft>>1)+1,dst ); + }else{ + int ncfft = nfft>>1; + int ncfft2 = nfft>>2; + Complex * rtw = real_twiddles(ncfft2); + + // use optimized mode for even real + fwd( dst, reinterpret_cast<const Complex*> (src), ncfft); + Complex dc = dst[0].real() + dst[0].imag(); + Complex nyquist = dst[0].real() - dst[0].imag(); + int k; + for ( k=1;k <= ncfft2 ; ++k ) { + Complex fpk = dst[k]; + Complex fpnk = conj(dst[ncfft-k]); + Complex f1k = fpk + fpnk; + Complex f2k = fpk - fpnk; + Complex tw= f2k * rtw[k-1]; + dst[k] = (f1k + tw) * Scalar(.5); + dst[ncfft-k] = conj(f1k -tw)*Scalar(.5); + } + dst[0] = dc; + dst[ncfft] = nyquist; } + } - inline - void fwd( Complex * dst,const Complex *src,int nfft) - { - get_plan(nfft,false).work(0, dst, src, 1,1); - } + // inverse complex-to-complex + inline + void inv(Complex * dst,const Complex *src,int nfft) + { + get_plan(nfft,true).work(0, dst, src, 1,1); + } - // real-to-complex forward FFT - // perform two FFTs of src even and src odd - // then twiddle to recombine them into the half-spectrum format - // then fill in the conjugate symmetric half - inline - void fwd( Complex * dst,const Scalar * src,int nfft) - { - if ( nfft&3 ) { - // use generic mode for odd - m_tmpBuf1.resize(nfft); - get_plan(nfft,false).work(0, &m_tmpBuf1[0], src, 1,1); - std::copy(m_tmpBuf1.begin(),m_tmpBuf1.begin()+(nfft>>1)+1,dst ); - }else{ - int ncfft = nfft>>1; - int ncfft2 = nfft>>2; - Complex * rtw = real_twiddles(ncfft2); - - // use optimized mode for even real - fwd( dst, reinterpret_cast<const Complex*> (src), ncfft); - Complex dc = dst[0].real() + dst[0].imag(); - Complex nyquist = dst[0].real() - dst[0].imag(); - int k; - for ( k=1;k <= ncfft2 ; ++k ) { - Complex fpk = dst[k]; - Complex fpnk = conj(dst[ncfft-k]); - Complex f1k = fpk + fpnk; - Complex f2k = fpk - fpnk; - Complex tw= f2k * rtw[k-1]; - dst[k] = (f1k + tw) * Scalar(.5); - dst[ncfft-k] = conj(f1k -tw)*Scalar(.5); - } - dst[0] = dc; - dst[ncfft] = nyquist; + // half-complex to scalar + inline + void inv( Scalar * dst,const Complex * src,int nfft) + { + if (nfft&3) { + m_tmpBuf1.resize(nfft); + m_tmpBuf2.resize(nfft); + std::copy(src,src+(nfft>>1)+1,m_tmpBuf1.begin() ); + for (int k=1;k<(nfft>>1)+1;++k) + m_tmpBuf1[nfft-k] = conj(m_tmpBuf1[k]); + inv(&m_tmpBuf2[0],&m_tmpBuf1[0],nfft); + for (int k=0;k<nfft;++k) + dst[k] = m_tmpBuf2[k].real(); + }else{ + // optimized version for multiple of 4 + int ncfft = nfft>>1; + int ncfft2 = nfft>>2; + Complex * rtw = real_twiddles(ncfft2); + m_tmpBuf1.resize(ncfft); + m_tmpBuf1[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() ); + for (int k = 1; k <= ncfft / 2; ++k) { + Complex fk = src[k]; + Complex fnkc = conj(src[ncfft-k]); + Complex fek = fk + fnkc; + Complex tmp = fk - fnkc; + Complex fok = tmp * conj(rtw[k-1]); + m_tmpBuf1[k] = fek + fok; + m_tmpBuf1[ncfft-k] = conj(fek - fok); } + get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_tmpBuf1[0], 1,1); } + } - // inverse complex-to-complex - inline - void inv(Complex * dst,const Complex *src,int nfft) - { - get_plan(nfft,true).work(0, dst, src, 1,1); - } + protected: + typedef ei_kiss_cpx_fft<Scalar> PlanData; + typedef std::map<int,PlanData> PlanMap; - // half-complex to scalar - inline - void inv( Scalar * dst,const Complex * src,int nfft) - { - if (nfft&3) { - m_tmpBuf1.resize(nfft); - m_tmpBuf2.resize(nfft); - std::copy(src,src+(nfft>>1)+1,m_tmpBuf1.begin() ); - for (int k=1;k<(nfft>>1)+1;++k) - m_tmpBuf1[nfft-k] = conj(m_tmpBuf1[k]); - inv(&m_tmpBuf2[0],&m_tmpBuf1[0],nfft); - for (int k=0;k<nfft;++k) - dst[k] = m_tmpBuf2[k].real(); - }else{ - // optimized version for multiple of 4 - int ncfft = nfft>>1; - int ncfft2 = nfft>>2; - Complex * rtw = real_twiddles(ncfft2); - m_tmpBuf1.resize(ncfft); - m_tmpBuf1[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() ); - for (int k = 1; k <= ncfft / 2; ++k) { - Complex fk = src[k]; - Complex fnkc = conj(src[ncfft-k]); - Complex fek = fk + fnkc; - Complex tmp = fk - fnkc; - Complex fok = tmp * conj(rtw[k-1]); - m_tmpBuf1[k] = fek + fok; - m_tmpBuf1[ncfft-k] = conj(fek - fok); - } - get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_tmpBuf1[0], 1,1); - } - } + PlanMap m_plans; + std::map<int, std::vector<Complex> > m_realTwiddles; + std::vector<Complex> m_tmpBuf1; + std::vector<Complex> m_tmpBuf2; - protected: - typedef ei_kiss_cpx_fft<Scalar> PlanData; - typedef std::map<int,PlanData> PlanMap; - - PlanMap m_plans; - std::map<int, std::vector<Complex> > m_realTwiddles; - std::vector<Complex> m_tmpBuf1; - std::vector<Complex> m_tmpBuf2; - - inline - int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; } - - inline - PlanData & get_plan(int nfft,bool inverse) - { - // TODO look for PlanKey(nfft, ! inverse) and conjugate the twiddles - PlanData & pd = m_plans[ PlanKey(nfft,inverse) ]; - if ( pd.m_twiddles.size() == 0 ) { - pd.make_twiddles(nfft,inverse); - pd.factorize(nfft); - } - return pd; + inline + int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; } + + inline + PlanData & get_plan(int nfft,bool inverse) + { + // TODO look for PlanKey(nfft, ! inverse) and conjugate the twiddles + PlanData & pd = m_plans[ PlanKey(nfft,inverse) ]; + if ( pd.m_twiddles.size() == 0 ) { + pd.make_twiddles(nfft,inverse); + pd.factorize(nfft); } + return pd; + } - inline - Complex * real_twiddles(int ncfft2) - { - std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there - if ( (int)twidref.size() != ncfft2 ) { - twidref.resize(ncfft2); - int ncfft= ncfft2<<1; - Scalar pi = acos( Scalar(-1) ); - for (int k=1;k<=ncfft2;++k) - twidref[k-1] = exp( Complex(0,-pi * ((double) (k) / ncfft + .5) ) ); - } - return &twidref[0]; + inline + Complex * real_twiddles(int ncfft2) + { + std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there + if ( (int)twidref.size() != ncfft2 ) { + twidref.resize(ncfft2); + int ncfft= ncfft2<<1; + Scalar pi = acos( Scalar(-1) ); + for (int k=1;k<=ncfft2;++k) + twidref[k-1] = exp( Complex(0,-pi * ((double) (k) / ncfft + .5) ) ); } - }; + return &twidref[0]; + } +}; + +/* vim: set filetype=cpp et sw=2 ts=2 ai: */ + |