%-------------------------------------------------------------------------------
%              ------< Problem 2.4 (Text, 5.4.2) >------
%-------------------------------------------------------------------------------
% Script File: hw2_4
% Plots the trigonometric interpolant of the Gauss Pallas data.
% Given that data, use CSInterp and CSeval to generate data at a grid with size
% some power of 2; then use that as the data to drive the DFT. Once this is 
% done, replace the call to DFT by a call to FFTRecur. Repeat for various powers
% and plot flop counts (and timings) for three cases (DFT, FFTRecur, CSInterp).
% Show a plot of the FFTRecur values against original data to establish
% correctness.
%-------------------------------------------------------------------------------
% Use some arbitrary data (here data from text, p. 102)
A = linspace(0,360,13)';
D=[ 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804 408]';
% now use CSInterp and CSeval to interpolate data to a grid that is 
% a power of 2; then use that as the original data, on which to run
% the three algorithms and perform benchmarking. Use powers of 2 in range
%  2^k: k = 3,...,8
%-------------------------------------------------------------------------------
disp('----------------------------------------------------------------')
disp('  k      n           DFT     FFTRecur     CSInterp       CSeval ')
disp('----------------------------------------------------------------')
for k = 3:1:6
   m=2^k;n=2*m;
   Avals = linspace(0,360,n+1)';
   size(Avals);
   F     = CSInterp(D(1:12));
   Fvals = CSeval(F,360,Avals(1:n));
%-------------------------------------------------------------------------------
%  flops(0);tic;
   y1 = DFT(Fvals,1);
%  f1 = flops;
   t1 = toc;
%  flops(0);
tic;
%  y2 = call to FFTRecur
%   f2 = flops;
   t2 = toc;
%   flops(0); 
tic;
   P3 = CSInterp(Fvals(1:n));
%  f3 = flops; 
   t3 = toc;
%-------------------------------------------------------------------------------
%  convert data from DFT, (FFTRecur) to data structures, as expected by
%  CSeval, using the formulas as given 
%  Note that the scaling by 1/m is performed inside DFT, but not inside 
%  FFTRecur, "to avoid unnessecary divisions by 2 at each iteration" 
%  The following comments marked by *** are not required! 
%  *** How could we scale inside FFTRecur? 
%  *** Why would that be more costly (although not by much!)?)
%  Note that if we run DFT with isign = -1 we get the Inverse Fourier Transform!
%  *** How would you change FFTRecur to produce the inverse transform?
   P1 = struct('a',real(y1(1:m+1)),'b',-imag(y1(2:m)));
   P1.a(1)=P1.a(1)/2;P1.a(m+1)=P1.a(m+1)/2;
%  P2 = struct('a',part of y2,'b',part of y2);
%-------------------------------------------------------------------------------
% Run CSeval on DFT (FFTRecur) generated coefficients to produce interpolated
% values for validation plot (do it each time for timing purposes, although
% only final values will be plotted)
%  flops(0);
tic;
   Bvals = linspace(0,360,360)';
   Gvals = CSeval(P1,360,Bvals);
%  Gvals = CSeval(P2,360,Bvals);
%   f4 = flops;
   t4 = toc;
disp(sprintf(' %2.0f   %4.0f    %10.0f   %10.0f   %10.0f   %10.0f',k,n,f1,f2,f3,f4))
end
disp('----------------------------------------------------------------')
%-------------------------------------------------------------------------------
plot(Bvals,Gvals,A,D,'o')
axis([-10 370 -200 1700])
set(gca,'xTick',linspace(0,360,13))
xlabel('Ascension (Degrees)')
ylabel('Declination (minutes)')
%-------------------------------------------------------------------------------
% OUTPUT of script: hw2_4
%-------------------------------------------------------------------------------
% >> hw2_4
% ----------------------------------------------------------------
%   k      n           DFT     FFTRecur     CSInterp       CSeval
% ----------------------------------------------------------------
%   3     16          3562         1005         5176        21975
%   4     32         13770         2469        31560        45023
%   5     64         54154         5845       213468        91119
%   6    128        214794        13493      1550874       183311
%   7    256        855562        30581     11797644       367695
%   8    512       3415050        68341     91904220       736463
% ----------------------------------------------------------------
% ----------------------------------------------------------------
% ----------------------------------------------------------------
% ----------------------------------------------------------------
%  functions called DFT, FFTRecur, CSeval, CSInterp
% ----------------------------------------------------------------
% ----------------------------------------------------------------
%   function y = DFT(x,isign)
% ----------------------------------------------------------------
% y = DFT(x,+/-1):  +1 input:point sample   -> output:freq. dom. coeffs.
%                   -1 input:freq.dom.coef. -> output:point values
% n = length(x);
% y = x(1)*ones(n,1);
% if n>1
%    m=n/2;
%    w=exp(isign*pi*sqrt(-1)/m);
%    v = w.^(0:n-1)';
%    for k = 2:n
%       z = rem((k-1)*(0:n-1)',n )+1;
%       y = y + v(z)*x(k);
%    end
%    if isign == 1
%       y = y/m;
%    end
% end
% ----------------------------------------------------------------
% ----------------------------------------------------------------
%  function y = FFTRecur(x)
% ----------------------------------------------------------------
% y = FFTRecur(x)
% y is the discrete Fourier transform of column vector x, n = 2^k
% 
% n = length(x);
% if n == 1
%    y = x;
% else
%    m = n/2;
%    w = exp(pi*sqrt(-1)/m);
%    yT = FFTRecur(x(1:2:n));
%    yB = FFTRecur(x(2:2:n));
%    d  = w.^(0:m-1)';
%    z  = d.*yB;
%    y  = [ yT+z ; yT-z ];
% end
% ----------------------------------------------------------------
% ----------------------------------------------------------------
%  function Fvals = CSeval(F,T,tvals)
% ----------------------------------------------------------------
% F is a structure: F('a',***,'b',***)
% F.a is a length m+1 column vector, F.b is a length m-1 column vector,
% T is a positive scalar, and tvals is a column vector.
% If 
%   F(t) = F.a(1) + F.a(2)*cos((2*pi/T)*t) +...+ F.a(m+1)*cos((2*m*pi/T)*t) +
%                   F.b(1)*sin((2*pi/T)*t) +...+ F.b(m-1)*cos((2*m*pi/T)*t)
%		   
% then Fvals = F(tvals).
% 
%  Fvals = zeros(length(tvals),1);
%  tau = (2*pi/T)*tvals;
%  for j = 0:length(F.a)-1, Fvals = Fvals + F.a(j+1)*cos(j*tau); end
%  for j = 1:length(F.b),   Fvals = Fvals + F.b(j  )*sin(j*tau); end
% ----------------------------------------------------------------
% ----------------------------------------------------------------
%  function F =CSInterp(f)
% ----------------------------------------------------------------
%  F =CSInterp(f)
% f 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
% tau = (pi/m)*(0:n-1)', then
%
%	f = F.a(1)*cos(0*tau) +...+ F.a(m+1)*cos( m   *tau) +
%	    F.b(1)*sin(  tau) +...+ F.b(m-1)*sin((m-1)*tau)
% 
% n = length(f);
% m = n/2;
% tau = (pi/m)*(0:n-1)';
% P = [];
% for j = 0:m,     P = [P cos(j*tau)]; end
% for j = 1:m-1,   P = [P sin(j*tau)]; end
% y = P\f;
% F = struct('a',y(1:m+1),'b',y(m+2:n));
% ----------------------------------------------------------------
