3
$\begingroup$

I want to use FFTs to do seriously accurate interpolation on band-limited data. To do that, I need to get a handle on the fundamental accuracy of the FFT() and IFFT() algorithms available. My idea is to look at IFFT(FFT(X)) which should be a unit transformation. I want to get the maximum difference between the input and output data down below the region of 1 part in 10^8 or even 10^9. Currently I'm achieving ~1 part in 10^6 using 64-bit float, but I don't know if that is in line with what is to be expected. I'm thinking that surely someone has done all this before. Does anybody know the answers? What can expect if I try using higher-precision data types? With 64-bit float speed is not a concern.

$\endgroup$
  • 1
    $\begingroup$ If you talk about a maximum error you need to specify the FFT length, because the error grows with length. Have you tried FFTW? It might be hard to get better than that. $\endgroup$ – Matt L. Dec 16 '15 at 21:31
  • $\begingroup$ Which FFT library or code are you using? $\endgroup$ – hotpaw2 Dec 16 '15 at 21:59
  • $\begingroup$ Mostly I've been using vDSP and a FFT length of 32768. $\endgroup$ – Richard M Dec 17 '15 at 15:13
1
$\begingroup$

Your mileage may vary, but in multiple environments on my x86_64 bit linux machine, I get much less error than 1 / 10^9.

julia:

julia> a=randn(32768); f=ifft(fft(a)); print(maximum(abs(f-a)))
1.782802105510316e-15

octave:

octave:1> a=randn(32768,1); f=ifft(fft(a)); disp(max(abs(f-a)))
1.8184e-15

numpy:

In [1]: %pylab
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [2]: a=randn(32768); f=fft.ifft(fft.fft(a)); print(max(abs(f-a)))
7.58800755928e-15

The default for all these environments is float64. If you want to share how the signal is bandlimited, that might shed more light on the problem.

UPDATE: To satisfy my curiosity, I tried this with a single sine-wav and basically get the same results:

julia> th=collect(0:32767)*(1000*pi/32767); a=sin(th); 
       f=ifft(fft(a)); print(maximum(abs(f-a)))
9.393934579042684e-16

Again, without more knowledge about the signal or the computation environment, this is hard to address.

$\endgroup$
  • 1
    $\begingroup$ Band limiting may indeed be a problem. Remember that the IFFT contains one bin per frequency band, and each bin holds 2x64 bits of data. (complex). A pure sine that fits the FFT period would be represented by just a single bin, 128 bits. Your random data is distributed over all bins (16385 x 128 bits). $\endgroup$ – MSalters Mar 1 '16 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.