Re: FFT module thoughts, was Re: [eigen] Eigen 3.0.5 Could NOT find FFTW

[ Thread Index | Date Index | More Archives ]

Good catch.  That explains why it worked properly when nfft was a power of 2.

On 04/05/2012 12:55 PM, Björn Piltz wrote:
The problem seems to be:

        scale(dst,Scalar(1./nfft),nfft); // scale the time series

which can be fixed as:

        scale(dst,Scalar(1.)/nfft,nfft); // scale the time series

That way you keep your precision.


2012/4/5 Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx>
Hi Mark.

Thank you for your prompt response.

I do not see problems with exp(), acos() or trig functions right now. Your code is exellent in taking care of generic Scalars.

For example, you use acos(Scalar(-1)), which works perfectly when Scalar=mpreal.
Same for exp() applied for std::complex<mpreal>, which is also defined and available to compiler.

Anyway, I'll double check correct function call tomorrow again.

I would really appreciate if you would provide some links on how to re-write make_twiddles using root finding approach.

Thank you for your time and efforts improving this world :).

On 05.04.2012, at 22:57, Mark Borgerding <mark@xxxxxxxxxxxxxx> wrote:

> I tried to keep this in mind when writing originally.  It does not handle it, but it might not be too hard.
> There are at least two Achille's heels,  both in kiss_cpx_fft::make_twiddles  ( file==ei_kissfft_impl.h )
> 1. The exp() function, as others have mentioned.
> 2. The calculation of pi, as done via acos(-1)
> Stepping back, what make_twiddles does is make one full cycle of a complex sinusoid with nfft points and unity magnitude.
> This problem easily decomposes into a primitive root-finding problem.  This can be done very efficiently in most cases for the FFT, since people tend to use a length that is a multiple of small primes (2,3,5)
> (Let me know if this is not clear and I can provide an example.)
> In other words, if we had a root-finding function
>  template<typename Scalar>
>  Scalar  nth_root( const RealScalar & x, int n); // maybe implemented via Newton-Raphson ?
> then make_twiddles could be specialized to use it for types other than complex<float|double|long double>
> The real_twiddles function also contains some trig calls -- these could derived in a similar way.
> On 04/05/2012 07:04 AM, Pavel Holoborodko wrote:
>> Yes, I have my own multiple precision complex exp() in std space, as
>> well as other math. functions. Seems to work fine in all other
>> contexts (including Eigen itself), but doesn't help with FFT :(
>> On Apr 5, 2012, at 19:52, "Björn Piltz"<bjornpiltz@xxxxxxxxxxxxxx>  wrote:
>>> I think std::exp() might be a candidate.
>>> Try adding
>>>     namespace std {  std::complex<Real>  exp(const std::complex<Real>&); }
>>> to your file to see if it's the culprit.
>>> Björn

Mail converted by MHonArc 2.6.19+