// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2009 Mark Borgerding mark a borgerding net // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
namespace Eigen {
namespace internal {
// This FFT implementation was derived from kissfft http:sourceforge.net/projects/kissfft // Copyright 2003-2009 Mark Borgerding
inlinevoid make_twiddles(int nfft, bool inverse)
{ using numext::sin; using numext::cos;
m_inverse = inverse;
m_twiddles.resize(nfft); double phinc = 0.25 * double(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)
{
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)
{ //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;
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;
}
}
/* 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 = static_cast<int>(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 += static_cast<int>(fstride) * k; if (twidx>=Norig) twidx-=Norig;
t=scratchbuf[q] * twiddles[twidx];
Fout[ k ] += t;
}
k += m;
}
}
}
};
// 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);
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
