function  [y energy] = fseries(t,func,M,T)
%  FSERIES
% general Fourier series evaluation
%
%      t     time values
%      func  function defining Fourier coefficients
%      M     number of terms to evaluate
%      T     period
%
%      If no output arguments, plot the time-signal reconstruction
%

if (nargin<4)
    T = 1.0;
end
if (nargin<3)
    M = 50;
end
f = 1/T;

% evaluate DC term, if NaN then sneak up on the value
c0 = func(0);
if (isnan(c0))
    c0 = func(1e-5);
end
energy = c0*conj(c0);
y = ones(size(t))*c0;
% now do the other terms in the series, pair up cosine and sine terms
for k=1:M
    arg = 2*pi*k*f*t;
    ck = func(k);
    ckm = func(-k);
    energy = energy + ck*conj(ck) + ckm*conj(ckm);
    y = y + (ck+ckm)*cos(arg)+ j*(ck-ckm)*sin(arg);
end
% calculate maximum absolute value of imaginary component
yi = imag(y);
tst = max(abs(yi));
% return the real component
y = real(y);
% do our own plotting
if (nargout<1)
    fprintf('energy (freq) %g (time) %g\n',energy, trapz(y.*y)/(length(t)-1));
    plot(t,y,'k','LineWidth',2);
    xlabel('time');
    ylabel('signal level');
    if (tst>1e-5)
        hold on;
        plot(t,yi);
        hold off;
    end
end

Published with MATLAB® 7.6