mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-27 07:29:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

Some improvements for kissfft from Martin Reinecke(pocketfft author):

Browse Source 
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.
This commit is contained in:
Guoqiang QI 2021-02-24 21:36:47 +00:00 committed by Rasmus Munk Larsen
parent a31effc3bc
commit f44197fabd
1 changed files with 39 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
unsupported/Eigen/src/FFT/ei_kissfft_impl.h
View File

@ -25,16 +25,47 @@ struct kiss_cpx_fft
std::vector<Complex> m_scratchBuf; std::vector<Complex> m_scratchBuf;
bool m_inverse; bool m_inverse;
inline inline void make_twiddles(int nfft, bool inverse)
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; Scalar c = Scalar(cos(i*8*phinc));
m_inverse = inverse; Scalar s = Scalar(sin(i*8*phinc));
m_twiddles.resize(nfft); m_twiddles[i] = Complex(c, s*flip);
Scalar phinc = (inverse?2:-2)* acos( (Scalar) -1) / nfft; m_twiddles[nfft-i] = Complex(c, -s*flip);
for (int i=0;i<nfft;++i)
m_twiddles[i] = exp( Complex(0,i*phinc) );
} }
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) void factorize(int nfft)
{ {