(* ::Package:: *)

function [a,b] = TrigInterpCVL(y)
% y is a column n vector and n = 2m.
% a is a column m+1 vector and b is a column m-1 vector so that if
%
% f(t) = (1/2)a(1)cos(0*t) + a(2)cos(t) + a(3)cos(2t) + ... + (1/2)a(m+1)cos(mt) +...
% b(1)sin(t) + b(2)sin(2t) + b(3)sin(3t) + ... + b(m-1)sin((m-1)t)
% then
% f(2*pi*i/n) = y(i+1), i=0:n-1
%
n = length(y);
m = n/2;
% We need z = (2/n)conj(Fn)*y for then a = z(0:m) and b = z(1:m-1)
% Take conjugates of both sides: conj(z) = (2/n)*Fn*conj(y) = (2/n)*Fn*y
conjz = (2/n)*fft(y);
a = real(conjz(1:m+1));
b = imag(conjz(2:m));

% -2 if did not play with conjugates properly.
% -3 if did not use fft