// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Mark Borgerding mark a borgerding net
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see .
// This FFT implementation was derived from kissfft http:sourceforge.net/projects/kissfft
// Copyright 2003-2009 Mark Borgerding
template
struct ei_kiss_cpx_fft
{
typedef _Scalar Scalar;
typedef std::complex Complex;
std::vector m_twiddles;
std::vector m_stageRadix;
std::vector m_stageRemainder;
std::vector m_scratchBuf;
bool m_inverse;
inline
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;in)
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
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;kreal() - Scalar(.5)*scratch[3].real() , Fout->imag() - Scalar(.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_twiddles.size());
Complex * scratchbuf = &m_scratchBuf[0];
for ( u=0; u(fstride) * k;
if (twidx>=Norig) twidx-=Norig;
t=scratchbuf[q] * twiddles[twidx];
Fout[ k ] += t;
}
k += m;
}
}
}
};
template
struct ei_kissfft_impl
{
typedef _Scalar Scalar;
typedef std::complex Complex;
void clear()
{
m_plans.clear();
m_realTwiddles.clear();
}
inline
void fwd( Complex * dst,const Complex *src,int nfft)
{
get_plan(nfft,false).work(0, dst, src, 1,1);
}
inline
void fwd2( Complex * dst,const Complex *src,int n0,int n1)
{
}
inline
void inv2( Complex * dst,const Complex *src,int n0,int n1)
{
}
// 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 (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;
}
}
// inverse complex-to-complex
inline
void inv(Complex * dst,const Complex *src,int nfft)
{
get_plan(nfft,true).work(0, dst, src, 1,1);
}
// 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>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(dst), &m_tmpBuf1[0], 1,1);
}
}
protected:
typedef ei_kiss_cpx_fft PlanData;
typedef std::map PlanMap;
PlanMap m_plans;
std::map > m_realTwiddles;
std::vector m_tmpBuf1;
std::vector m_tmpBuf2;
inline
int PlanKey(int nfft, bool isinverse) const { return (nfft<<1) | int(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 & 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 * (Scalar(k) / ncfft + Scalar(.5)) ) );
}
return &twidref[0];
}
};
/* vim: set filetype=cpp et sw=2 ts=2 ai: */