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

On Fri, Apr 6, 2012 at 1:55 AM, Björn Piltz <bjornpiltz@xxxxxxxxxxxxxx> 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.

Björn

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
>>>
>>
>
>
>