function [A,f] = spectrum(t,x,fmax,Tstart,Tstop);
% function [A] = spectrum(t,x,fmax,Tstart,Tstop);
%
% Calculates discrete fourier transform X for signal x(t)
% Extracts Fourier coefficients a and b from transform X 
% Converts to A(k)*cos(w*t + arg(k)), where k = 0,1,2... fmax/50.
%
% Note: Signal don't need to be sampled at equidistant steps
%
% Tstart and Tstop in milliseconds [ms]

dt = 1e-6;                                                 % [s]

[xi,ti] = interpolate(x,t,dt,Tstart/1000,Tstop/1000);

N = length(xi);

X = fft(xi);       % Get Fourier transform

a = 2*real(X)/N;   % Coeffs before cos(2*pi*k*x(n)/(N*dt))
b = 2*imag(X)/N;   % Coeffs before sine
Atemp = sqrt(a.^2 + b.^2);

A = Atemp(1:fmax/20+1);
f = [0:20:fmax]';

bar(f,A);
grid;
axis([min(f) max(f) min(A) max(A)]);