"alex singh" <singhalexander@n_o_s_p_a_m.gmail.com> wrote in message 
news:gJudnZ3S7teqR1bSnZ2dnUVZ_tadnZ2d@giganews.com...
> Hi, I was hoping someone could help me understand the effects of the
> numerical precision in Matlab's FFT. I setup a very simple code that
> perform the FFT and IFFT and compare my results to an FFT follwoed by a
> "manual" IFFT, I see a very strange increasing behavioir in the error and
> was wondering what this was due to. The fact that it increases with the
> sample number leads me to believe it's related to the index but I am not
> sure how I can mitigate this problem.

The trend you see is probably due to the build up in error when you 
calculate sin/cos in:
    yr = yr + xk(k)*exp(1i*2*pi*[0:N-1].'*(k-1)/N);
Matlab will try to calculate sin/cos for arguments with a magnitude of ~2^26 
* 2 * pi.

Try for example:
 cos(0.5) - cos(0.5 + N^2*2*pi)
It should be zero but it's not

/rgds 

kevin <kevinjmcgee@netscape.net> wrote:
> oops - I meant N log (base 2) N in the above, not 'ln'.

Big-Oh notation ignores constant factors, so any log base
will work. That is why it is usually just written log.

Is PL/I the only programming language with a log2() function?

-- glen

oops - I meant N log (base 2) N in the above, not 'ln'.

Kevin

On Jun 3, 5:42&#4294967295;pm, "alex singh" <singhalexander@n_o_s_p_a_m.gmail.com>
wrote:
> Hi, I was hoping someone could help me understand the effects of the
> numerical precision in Matlab's FFT. I setup a very simple code that
> perform the FFT and IFFT and compare my results to an FFT follwoed by a
> "manual" IFFT, I see a very strange increasing behavioir in the error and
> was wondering what this was due to. The fact that it increases with the
> sample number leads me to believe it's related to the index but I am not
> sure how I can mitigate this problem.
>
> CODE:
>
> N = 2^13;
> x = rand(N,1);
> xk = fft(x,N);
> y = ifft(xk,N);
> yr = zeros(N,1);
> for k=1:N
> &#4294967295; &#4294967295;yr = yr + xk(k)*exp(1i*2*pi*[0:N-1].'*(k-1)/N);
> end
> figure;
> plot(abs(yr/N-x))
> title 'hand ifft'
> figure;
> plot(abs(x-y))
> title 'ifft'
> Could it be the constants values? how does fftw actually handle these
> numbers.
>
> Thanks for the help.
>
> Alex


Well, I don't do MATLAB, so I can't help much with your code.  And
it's late here, so I won't write much.  But MATLAB FFT's are based on
FFTW, and the FFTW web site might tell you some of what you need to
know:

http://www.fftw.org/accuracy/comments.html

FFT error analysis is a well studied field going back to the 1960's.
There are literally several dozen research papers that deal with it
(for fixed, block floating and floating point) .  Much depends on how
the algorithm is coded, the accuracy of the twiddles used, etc.

No surprise that the DFT versus FFT would show different errors - the
DFT version should be worse than the FFT.  That's because the FFT gets
its results using fewer computations than a DFT (N ln N versus N
squared).

I'm not sure what you mean when you say 'related to the index.'  The
error should depend on N, of course, but FFT errors grow very slowly
relative to it (I believe the relationship is described in the above
link)

Kevin McGee

Hi, I was hoping someone could help me understand the effects of the
numerical precision in Matlab's FFT. I setup a very simple code that
perform the FFT and IFFT and compare my results to an FFT follwoed by a
"manual" IFFT, I see a very strange increasing behavioir in the error and
was wondering what this was due to. The fact that it increases with the
sample number leads me to believe it's related to the index but I am not
sure how I can mitigate this problem.

CODE:

N = 2^13;
x = rand(N,1);
xk = fft(x,N);
y = ifft(xk,N);
yr = zeros(N,1);
for k=1:N
   yr = yr + xk(k)*exp(1i*2*pi*[0:N-1].'*(k-1)/N);
end
figure;
plot(abs(yr/N-x))
title 'hand ifft'
figure;
plot(abs(x-y))
title 'ifft'
Could it be the constants values? how does fftw actually handle these
numbers.

Thanks for the help.

Alex