moved real-half-spectrum reflection into Eigen::FFT

This commit is contained in:
Mark Borgerding 2009-10-30 20:26:30 -04:00
parent a26b729cc9
commit d659fd9b14
3 changed files with 22 additions and 15 deletions

View File

@ -88,6 +88,8 @@ class FFT
void fwd( Complex * dst, const Scalar * src, int nfft) void fwd( Complex * dst, const Scalar * src, int nfft)
{ {
m_impl.fwd(dst,src,nfft); m_impl.fwd(dst,src,nfft);
if ( HasFlag(HalfSpectrum) == false)
ReflectSpectrum(dst,nfft);
} }
void fwd( Complex * dst, const Complex * src, int nfft) void fwd( Complex * dst, const Complex * src, int nfft)
@ -115,7 +117,11 @@ class FFT
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
dst.derived().resize( src.size() );
if ( NumTraits< typename InputDerived::Scalar >::IsComplex == 0 && HasFlag(HalfSpectrum) )
dst.derived().resize( (src.size()>>1)+1);
else
dst.derived().resize(src.size());
fwd( &dst[0],&src[0],src.size() ); fwd( &dst[0],&src[0],src.size() );
} }
@ -160,7 +166,6 @@ class FFT
inv( &dst[0],&src[0],dst.size() ); inv( &dst[0],&src[0],dst.size() );
} }
// TODO: multi-dimensional FFTs // TODO: multi-dimensional FFTs
// TODO: handle Eigen MatrixBase // TODO: handle Eigen MatrixBase
@ -176,6 +181,13 @@ class FFT
*x++ *= s; *x++ *= s;
} }
void ReflectSpectrum(Complex * freq,int nfft)
{
int nhbins=(nfft>>1)+1;
for (int k=nhbins;k < nfft; ++k )
freq[k] = conj(freq[nfft-k]);
}
impl_type m_impl; impl_type m_impl;
int m_flag; int m_flag;
}; };

View File

@ -177,9 +177,6 @@
void fwd( Complex * dst,const Scalar * src,int nfft) void fwd( Complex * dst,const Scalar * src,int nfft)
{ {
get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft); get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft);
int nhbins=(nfft>>1)+1;
for (int k=nhbins;k < nfft; ++k )
dst[k] = conj(dst[nfft-k]);
} }
// inverse complex-to-complex // inverse complex-to-complex

View File

@ -322,8 +322,6 @@
// place conjugate-symmetric half at the end for completeness // place conjugate-symmetric half at the end for completeness
// TODO: make this configurable ( opt-out ) // TODO: make this configurable ( opt-out )
for ( k=1;k < ncfft ; ++k )
dst[nfft-k] = conj(dst[k]);
dst[0] = dc; dst[0] = dc;
dst[ncfft] = nyquist; dst[ncfft] = nyquist;
} }