function F = CSinterpFast(fVals)
% fVals is a column n vector and n = 2m.
% F.a is a column m+1 vector and F.b is a column m-1 vector so that if
%
% f(t) = F.a(1)cos(0*t) + F.a(2)cos(t) + F.a(3)cos(2t) + ... + F.a(m+1)cos(mt) +...
% F.b(1)sin(t) + F.b(2)sin(2t) + F.b(3)sin(3t) + ... + F.b(m-1)sin((m-1)t)
% then
% f(2*pi*i/n) = fVals(i+1), i=0:n-1
%
n = length(fVals);
m = n/2;

% Solve Py = fVals where P = CSmatrix(n)
% I.e., solve (P'P)y = P'*fVals. Use the fact that D = P'*P is diagonal.
% In particular D = m*eye except D(1,1) = D(m+1,m+1) = n

% Form b = P'*fVals
b = zeros(n,1);
t = (0:n-1)'*(2*pi/n);
for j=0:m
b(j+1) = cos(j*t)'*fVals;
end
for j=1:m-1
b(j+m+1) = sin(j*t)'*fVals;
end
% Now solve Dy = b...
y = b/m; y(1) = y(1)/2; y(m+1) = y(m+1)/2;
% Package the coefficients...
F = struct('a',y(1:m+1),'b',y(m+2:n));