// Fourier Integral of f, here A(w) = 0, B(w) = (sin(w) + w cos(w))/w^2
function y = fi(x, n)
// Avoid w = 0
w = 0.001:0.01:n;
m = size(x,2);
y = x;
for i = 1:m, // for each x, compute the integral with resp to w
	fw = (sin(w) - w .* cos(w)) ./ (w .^2)  .*  sin( w * x(i));
	y(i) = inttrap(w,fw);
end; // ends the for statement
endfunction; // ends the function

t=-3:0.01:3;

// The function f is pi x/2 for -1 < x < 1
f=(%pi/2)*t;
f(find(t<-1)) = 0;
f(find(t>1)) = 0;

scf(0);
clf();
f1=fi(t,%pi);
f2=fi(t,2*%pi);
f4=fi(t,4*%pi);
f10=fi(t,10*%pi);

plot(t,f,'k-');
plot(t,f1,'m-');
plot(t,f2,'b-');
plot(t,f4,'g-');
plot(t,f10,'r-');
xtitle('f in black, N = Pi in magneta, N = 2Pi in blue, N = 4Pi in green, N = 10Pi in red','x','f');

// draw the graph in a postscript file
xs2eps(0,'fi.eps');