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?

Solution 1

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.

alt text http://files.droplr.com/files/35740123/XUl90.fft.png

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))

Solution 2

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!

Comments

[Clément]

Olivier, you rock :) That's exactly what I needed, Thanks!

[Clément]

Great, thanks! I got it now. Do you know of good documentation sources about fft?

[Clément]

Although, what does the use of the negative -cutoff subscript do?

[Olivier Verdier]

CFP, glad you like it! -cutoff means the "cutoff to last" index, i.e., -1 means the last index. So the slice [size/2,-cutoff] means leave everything from half, except the cutoff last. A neater way would have been: fourier[cutoff+1:-cutoff]=0.

[Martin Thoma]

I had to use scipy.fftpack instead of scipy.signal. And I got /usr/local/lib/python2.7/dist-packages/numpy/core/numeric.py‌​:460: ComplexWarning: Casting complex values to real discards the imaginary part return array(a, dtype, copy=False, order=order). And it did not plot anything (you should probably add show()).

[Martin Thoma]

Why do you calculate size?