function [alfa,beta,gamma,Predict] = TrigLS(tau,f,P,n)
% tau and f are column m-vectors, P is a positive real, and
% n is a positive integer (2n+1<=m).
% Determines alfa (scalar), beta (column n-vector), and gamma (column n-vector)
% so that if
%
% F(t) = alfa + sum_{k=1:n} of (beta(k)cos(2pi*k*t/P) + gamma(k)*sin(2*pi*k*t/P) )
% then
% sum_{i=1:m} of ( F(tau(i)) - f(i) )^2
%
% is minimized. In other words, the data (tau(1),f(1)),...,(tau(m),f(m)) is being
% fit in the least squares sense with a sum of sines and cosines.
% Predict is the column m-vector F(tau), i.e., Predict(i) = F(tau(i)),
% i=1:m.

% Set up an LS problem min norm(A*x-b,2) where x(1) = alfa,
% x(2:n+1) = beta, and x(n+2:2n+1) = gamma...
m = length(f);
b = f;
A = zeros(m,2*n+1);
A(:,1) = ones(m,1);
for k=1:n
tVals = (2*pi*k/P)*tau;
A(:,k+1) = cos(tVals);
A(:,n+k+1) = sin(tVals);
end
% Or...
% A = [ones(m,1) cos((2*pi/P)*tau*(1:n)) sin((2*pi/P)*tau*(1:n))];

% Now solve and take apart...
x = A\b;
alfa = x(1);
beta = x(2:n+1);
gamma = x(n+2:2*n+1);
Predict = A*x;