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.
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> a=randn(32768); f=ifft(fft(a)); print(maximum(abs(f-a))) 1.782802105510316e-15
octave:1> a=randn(32768,1); f=ifft(fft(a)); disp(max(abs(f-a))) 1.8184e-15
In : %pylab Using matplotlib backend: Qt4Agg Populating the interactive namespace from numpy and matplotlib In : 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.