% script: DFTtest.m
% description: perform some approximation experiments with complex
% Fourier polynomials on the function myfcn defined externally.

% programmer: t. shores
% last rev.: 2/14/16
 
% Octave and Matlab compatible file
n = 4 % highest frequency in Fourier approximation
N = 2*n + 1 % account for Nyquist sampling rate
a = zeros(1,n+1); 
b = zeros(1,n);
c = zeros(1,N);
% construct the coefficients for myfcn 
xnodes = linspace(0,2*pi,N+1); 
xnodes = xnodes(1:N);  % trim off 2*pi
f = myfcn(xnodes);
% dispose of 0th coefficient
a(n+1) = 2*sum(f)/N;
% now the rest of real Fourier coefs
for k = 1:n 
  a(k) = 2*f*cos(k*xnodes)'/N; 
  b(k) = 2*f*sin(k*xnodes)'/N;
end
% next construct the complex coefficients
% with a suitable DFT matrix
zeta = exp(-2*pi*i/N);
z = zeta.^(-n:n);
F = fliplr(vander(z))/N;
c = F*f.';
% ok, let's compare
ac = zeros(1,n+1);
% a(n+1) is exceptional Fourier a_0 of index 0
ac(n+1) = 2*c(n+1);
ac(1:n) = c(n+2:N) + c(n:-1:1)
bc = i*(c(n+2:N) - c(n:-1:1) ).';
disp('| a_0 - 2*c_0 | :') ,disp(abs(a(n+1) - ac(n+1)))
% the a(k) are Fourier a_k, k>0
disp('|| a_k - (c_k + c_{-k}) || :'), disp(norm(a(1:n)-ac(1:n),inf))
% the b(k) are Fourier b_k, k>0
disp(' || b_k - (c_k - c_{-k}) || :'), disp(norm(b - bc,inf))
