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!

Answer 2

You are probably doing it wrong. You remove all the "negative" frequencies in your code. You should keep both positive and negative low frequencies. Here is a code in python and the result. The plot has the right scale.

The code:

from __future__ import division

from scipy.signal import fft, ifft
import numpy as np

def approximate(signal, cutoff):
    fourier = fft(signal)
    size = len(signal)
    # remove all frequencies except ground + offset positive, and offset negative:
    fourier[1+cutoff:-cutoff] = 0
    return ifft(fourier)

def quad(x):
    return x**2

from pylab import plot

X = np.arange(-4096,4096)/64
Y = quad(X)

plot(X,Y)
plot(X,approximate(Y,3))