diff options
| author |  Guoqiang QI <425418567@qq.com> | 2021-02-24 21:36:47 +0000 |
|---|---|---|
| committer |  Rasmus Munk Larsen <rmlarsen@google.com> | 2021-02-24 21:36:47 +0000 |
| commit | f44197fabdfeb6008fba543e0a12959dc20735ae (patch) | |
| tree | 0520a92d63af484b8fef04395ef50b8c0679a230 /unsupported | |
| parent | a31effc3bca2f0924752caeebfd6f61f7edf9a43 (diff) | |
Some improvements for kissfft from Martin Reinecke(pocketfft author):
1.Only computing about half of the factors and use complex conjugate symmetry for the rest instead of all to save time.
2.All twiddles are calculated in double because that gives the maximum achievable precision when doing float transforms.
3.Reducing all angles to the range 0<angle<pi/4 which gives even more precision.
Diffstat (limited to 'unsupported')
| -rw-r--r-- | unsupported/Eigen/src/FFT/ei_kissfft_impl.h | 47 |
1 files changed, 39 insertions, 8 deletions
diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h index 9bdc012dd..303f85f13 100644 --- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h +++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h @@ -25,16 +25,47 @@ struct kiss_cpx_fft std::vector<Complex> m_scratchBuf; bool m_inverse; - inline - void make_twiddles(int nfft,bool inverse) + inline void make_twiddles(int nfft, bool inverse) + { + using numext::sin; + using numext::cos; + m_inverse = inverse; + m_twiddles.resize(nfft); + double phinc = 0.25 * EIGEN_PI / nfft; + Scalar flip = inverse ? Scalar(1) : Scalar(-1); + m_twiddles[0] = Complex(Scalar(1), Scalar(0)); + if ((nfft&1)==0) + m_twiddles[nfft/2] = Complex(Scalar(-1), Scalar(0)); + int i=1; + for (;i*8<nfft;++i) { - using std::acos; - 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) ); + Scalar c = Scalar(cos(i*8*phinc)); + Scalar s = Scalar(sin(i*8*phinc)); + m_twiddles[i] = Complex(c, s*flip); + m_twiddles[nfft-i] = Complex(c, -s*flip); + } + for (;i*4<nfft;++i) + { + Scalar c = Scalar(cos((2*nfft-8*i)*phinc)); + Scalar s = Scalar(sin((2*nfft-8*i)*phinc)); + m_twiddles[i] = Complex(s, c*flip); + m_twiddles[nfft-i] = Complex(s, -c*flip); } + for (;i*8<3*nfft;++i) + { + Scalar c = Scalar(cos((8*i-2*nfft)*phinc)); + Scalar s = Scalar(sin((8*i-2*nfft)*phinc)); + m_twiddles[i] = Complex(-s, c*flip); + m_twiddles[nfft-i] = Complex(-s, -c*flip); + } + for (;i*2<nfft;++i) + { + Scalar c = Scalar(cos((4*nfft-8*i)*phinc)); + Scalar s = Scalar(sin((4*nfft-8*i)*phinc)); + m_twiddles[i] = Complex(-c, s*flip); + m_twiddles[nfft-i] = Complex(-c, -s*flip); + } + } void factorize(int nfft) { |