function [c] = fseries( fun, a, b, N)
%FSERIES Fourier series

%	Giuseppe Rodriguez,  15-12-14
%	Email: rodriguez@unica.it

T = b-a;
omega = 2*pi/T;

t = [0:N-1]'/N * (b-a) + a;
f = fun(t);
fs = fftshift(f);
c = 1/N * fft(fs);
c = fftshift(c);

