Using GNU Octave FFT functions

Question

I'm playing with octave's fft functions, and I can't really figure out how to scale their output: I use the following (very short) code to approximate a function:

function y = f(x)
    y = x .^ 2;
endfunction;

X=[-4096:4095]/64;
Y = f(X);
# plot(X, Y);

F = fft(Y);
S = [0:2047]/2048;

function points = approximate(input, count)
    size    = size(input)(2);
    fourier = [fft(input)(1:count) zeros(1, size-count)];
    points  = ifft(fourier);
endfunction;

Y = f(X); plot(X, Y, X, approximate(Y, 10));

Basically, what it does is take a function, compute the image of an interval, fft-it, then keep a few harmonics, and ifft the result. Yet I get a plot that is vertically compressed (the vertical scale of the output is wrong). Any ideas?

Answer 1

You are throwing out the second half of the transform. The transform is Hermitian symmetric for real-valued inputs and you have to keep those lines. Try this:

function points = approximate(inp, count)
    fourier = fft(inp);
    fourier((count+1):(length(fourier)-count+1)) = 0;
    points  = real(ifft(fourier)); %# max(imag(ifft(fourier))) should be around eps(real(...))
endfunction;

The inverse transform will invariably have some tiny imaginary part due to numerical computation error, hence the real extraction.

Note that input and size are keywords in Octave; clobbering them with your own variables is a good way to get really weird bugs down the road!