% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fourierSC
% ____________________________________________________________________________
% 18.385 (Nonlinear Dynamical Systems).         MIT, Fall 1999. R. R. Rosales.
%
% Demo: To illustrate the convergence of the Fourier Series of a 2*pi periodic
% function F = F(x). The Fourier Series here is defined by
%     
%  F(x) = cF_0 + Sum_{n=1:infinity} [ cF_n * cos(n*x) + sF_n * sin(n*x) ]. (1)
%
% The demo computes the Fourier Coefficients cF_n and sF_n  up to some n = Np,
% that is specified by the user via screen input.   The demo has several built
% in examples of functions F(x) --- which the user can enter as an OPTION, via
% screen input. Arbitrary periodic functions F = F(x) are allowed,  provided a
% file FSFun.m (with a function script) is present in the directory from which
% the user is running MatLab. The MatLab statement y = FSFun(x)  should return
% the values of F(x) in the array y [i.e. y(n)=F(x(n))],  for any array x with
% x(n) in [0, 2*pi].       The function in FSFun will be used if OPTION = 0 is
% picked when running the demo (a default is provided).
%                               ^^^^^^^^^^^^^^^^^^^^^
% A simple example of a FSFun.m file is:
%                                       function y = FSFun(x)
%                                       y = exp(-10*(x-pi).^2);
% Note the use of the vector operand .^, needed to get a vector valued answer.
%
% >>>> Check help FSFun for important details about how to define the function
% when dealing with discontinuities!
%
% >>>> NOTE: help FSFun will NOT work if the current directory for MatLab has
% a user constructed  FSFun.m file already in it.
% ____________________________________________________________________________
%
% The demo does the following:
%   a) Computes the Fourier Coefficients and stores them in the 1 by Np arrays
%      sF and cF. The mean value cF_0 is stored in the variable cF_0.
%   b) Plots the sine and cosine Fourier coefficients as functions of n.
%   c) Shows how F is approximated as terms are added to the series. Plots are
%      shown of the approximations to (1) above given by:
%
%               cF_0 + Sum_{n=1:N} [ cF_n * cos(n*x) + sF_n * sin(n*x) ],  (2)
%
%      for N=0:Nskip:Np  (including Np, even if Np is not a multiple of Nskip)
%      and compared with the exact function. For N=0 the plot is just the mean
%      value of the function F (i.e. cF_0). The values of Nskip and Np must be
%      picked by the user (via requested screen input).
%   d) Plots the error in the approximation above as a function of N. Error is
%      in the array goof,  where goof(n) = error after N = n-1 coefficients in
%      the sum, with N=0:Np.  The error is measured pointwise in the numerical
%      grid x and then normalized by the range of F(x) = max(F) - min(F). This
%      except near discontinuities in F(x). Near the discontinuities of F(x) a
%      refined computational grid is introduced, and a calculation of how well
%      the Fourier series approximation does the jump is done.  The error here
%      is defined simply as the difference between the two limits of F, at the
%      discontinuity, and the max and min achieved by the approximation there,
%      again: normalized by the range of F.
%   e) Plots the power spectrum as a function of n,   using semilog and loglog
%      scales. Here the power spectrum is defined by
%   
%           power(n+1) = sqrt((cF_n)^2 + (sF_n)^2)   for n=1:Np            (3)
%           power(1)   = abs(cF_0)
%
%      The decay rate of the power spectrum as n grows is directly related to,
%      for example, how well the Fourier Series converges.
% ____________________________________________________________________________
%
% The following variables are left in memory -- corresponding to the last F(x)
% done -- at the end of a run.  NOTE: if there are no discontinuities in F(x),
% then some of the variables that have to do with discontinuities will not  be
% defined. Specifically: xRef, recRef and nwidth;  while Idis and Jump will be
% null arrays.
%
% Np     number of Fourier coefficients computed.
% Nskip  approximations to the function F are plotted with N terms in the sums
%        --- with N = 0, Nskip, 2*Nskip, .... k*Nskip, Np.  Because of the way
%        the code works, it is a good idea to avoid taking Nskip too large ---
%        say, not much bigger than 50. Memory usage goes up with Nskip and the
%        code performance deteriorates.
% dx     grid spacing (defined in the code).
% x      row array with space mesh. Namely: x=dx:dx:2*pi, for some dx.
% fct    row array with the values  F(x).
% rec    row array with the Fourier Series approximation to F(x) --- as in (2)
%        above --- with N = Np.
% xRef   refined grid with points near discontinuities in F(x).  This grid has
%        a spacing of dx/nRef --- as opposed to dx for x --- and extends  (for
%        each discontinuity point x_dc)   a distance pl*dx to the left of x_dc
%        and a distance pr*dx to the right of x_dc.    Thus each discontinuity
%        has nwidth = nRef*(pl+pr)+1 associated with it in xRef.
% recRef row array with values of the Fourier Series approximation at xRef.
%        NOTE: the plots done by the script DO NOT include the data in recRef,
%              which is used only for the error calculation.
% pl     parameter defining xRef (defined in the code).
% pr     """"""""""""""""""""""""""""""""""""""""""""""
% nwidth """"""""""""""""""""""""""""""""""""""""""""""
% nRef   """"""""""""""""""""""""""""""""""""""""""""""
% Jump   array with the indexes j in x=x(j) that correspond to points in xRef.
% Idis   indexes in grid x where discontinuities (approximately) occur.  That
%        is: the discontinuities are at or near x(Idis(n)).
% DtR    trigger value for discontinuities. A discontinuity at x(j) is assumed
%        when abs(fct(j) - fct(j-1)) > DtR*dx a
% power  row array (size Np+1), with the power spectrum. See (3) above. 
% goof   row array (size Np+1), with the errors in the Fourier Sums --- in (2)
%        above --- as function of N.   The errors are relative to the range of
%        F(x), defined by range = max(F) - min(F).
% cF_0   mean value of F.
% cF     row array with cosine Fourier coefficients for n=1:Np.
% sF     row array with sine Fourier coefficients for n=1:Np.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
VERSION = 2;
fprintf('\n There are two versions: ')
fprintf('\n            fancy (with buttons);VERSION = 1,')
fprintf('\n         and pedestrian;         VERSION = 2. \n')
VERSION = input('                           ENTER VERSION = ');
if VERSION == 1
   VERSION = 1;
else
   VERSION = 2;
end
%========== Vector sizes and useful text variables: ==========================
M    = 2048;
N    = 2*M;         % Size of vectors in space. In the code dx = 2*pi/N.
DtR  = 75;          % Trigger: define discontinuities where dF/dx > DtR
pl   = 31;          % Buffer size to left  of a discontinuity is pl*dx
pr   = 30;          % Buffer size to right of a discontinuity is pr*dx
nRef = 3;           % Refined grid near discontinuity has size dx/nRef
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Notes here on the choice of pl, pr and nRef:
%
% A) Near a discontinuity we want to exclude a region of (at least) 1/8 of the
%    shortest wave-length in the sum,   so we do not end up computing as error
%    the zone where the Fourier Series connects: the mean value of the jump at
%    at the discontinuity with the two limit values of the function.
%
%    When we sum Np terms, the shortest wave-length is 2*pi/Np,   thus we want
%    pl*dx to be at least 1/8 of this. That is: pl = N/8*Np. Thus, for pl = 30
%    and N = 4096, things will work fine for Np > 17. A similar argument works
%    for pr.
%
%    Thus, with the numbers N=4096 and pl~pr~30,  error calculations should be
%    fine for sums with more than (say) 20 terms.     Below this it's not even
%    worth worrying,     as the errors at discontinuities will be fairly large
%    anyway.
%
% B) The second point is how much do we need to refine. If we want the refined
%    grid to have at least 25 points per period at the maximum number of terms
%    to be added (Np = 512), we need 25*dx/nRef = 2*pi/512,       which yields
%    nRef = 25*512/N = 3.125 for N = 4096.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FN   = 'FontName';
HNB  = 'HelveticaNarrowBold';
FS   = 'FontSize';
LW   = 'LineWidth';
MS   = 'MarkerSize';
%
%========== Enter inputs needed and print on-screen info: ====================
%
fprintf('\n The exact function F(x) and the Fourier approximations: \n')
fprintf('\n      cF_0 + Sum_{n=1:N}[cF_n*cos(n*x) + sF_n*sin(n*x)], \n')
fprintf('\n with N = 0:Nskip:Np will be plotted (N=0 is the mean of F).\n')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
myFLAG = 1;         % To allow continuous execution. ////////////////////////#
while myFLAG > 0.5  %
%
Np    = input('               Enter Np, with 1 < Np < 512 .......... Np = ');
Nskip = input('               Enter Nskip (with Nskip < 60) ..... Nskip = ');
Np = floor(Np); Nskip = floor(Nskip);
Np = min(max([Np, 1]), 512); 
Nskip = min(max([Nskip, 1]), 60);
%
%========== Construct array of x's and calculate function ====================
%
dx = 2*pi/N;
x  = 0:dx:2*pi-dx;
fprintf('\n NOW CHOOSE THE FUNCTION F.\n')
if VERSION == 1
   FSoption
   elseif VERSION == 2
   FSoptionP
end
%
%========== Now detect the indexes where discontinuities occur, etc. =========
%
% We need to do this to have a decent error evaluation  at the points  where F
% is discontinuous, because the Fourier series oscillates there and point-wise
% error evaluation leads to misleading results.
%
% We ASSUME that discontinuities occupy at most 3 points in the grid (that is:
% "large" jumps in fct can occur in no more than two  consecutive intervals in
% the index j.    In fact, we ASSUME that the function whose Fourier series is
% being computed is piecewise continuous  and that the evaluation near a point
% of discontinuity (say, at x = x_d) is done by:
%
%  (A) F(x) defined by continuous expressions for both x > x_d and x < x_d;
%      with some value assigned at x = x_d.
%  (B) The value assigned at x_d satisfies:
%             F(x_d) = Fl  or  F(x_d) = Fr or F(x_d) = 0.5*(Fl + Fr),
%      where: Fl and Fr are the limits of F(x) from the left and right at x_d.
%      When F(x_d) = Fl or F(x_d) = Fr,  the "large" jump in fct will occur in
%      just one grid interval (from some x_j to x_(j+1)).   Otherwise the jump
%      may end up spread over two consecutive grid intervals --- note that the
%      periodicity condition identifies N+1 with 1.
%
% We detect discontinuities by checking for large "jumps" in fct,  where large
% means that the numerical derivative is larger than some threshold DtR.  Each
% discontinuity is labeled by one index: the right one in the jump if two grid
% points are involved and the middle one in the three grid points case.  These
% indexes are in the array Idis.
%
Temp = [abs(fct(1)-fct(N)), abs(fct(2:N)-fct(1:N-1))]/dx;
Idis = find((Temp > DtR) & (Temp([N, 1:N-1]) <= DtR));
Ndis = length(Idis);
%
% At each discontinuity:  a region of width pl*dx to the left and pr*dx to the
% right will be done carefully.   A refined grid will be constructed there for
% this, in the array xRef. The indexes for the points in x affected are in the
% array Jump. The arrays Maxdisf and Mindisf contain the approximate values of
% the limits for the function on each discontinuity.    The refining factor in
% xRef is nRef, i.e. dx ---> dx/nRef.
%
Jump = [];
if Ndis >= 1
   xRef = [];
   % ---------------------- Column j in indx is [(j+pp):(-1):(j-qq)] mod N
   pp = 1; qq = 1;
   indx = toeplitz([pp+1:-1:1, N:-1:N-qq+1], [pp+1:N, 1:pp]);
   for nn=1:Ndis
       Maxdisf(nn) = max(fct(indx(:, Idis(nn))));
       Mindisf(nn) = min(fct(indx(:, Idis(nn))));
   end
   % ---------------------- Column j in indx is [(j+pp):(-1):(j-qq)] mod N
   pp = pr; qq = pl;
   indx = toeplitz([pp+1:-1:1, N:-1:N-qq+1], [pp+1:N, 1:pp]);
   for nn=1:Ndis
       Jump = [Jump, indx(:, Idis(nn))'];
       xRef = [xRef, (x(Idis(nn))-pl*dx):(dx/nRef):(x(Idis(nn))+pr*dx)];
   end
   nwidth = nRef*(pl+pr)+1;  % Number of points in xRef at each discontinuity.
   Jump = sort(Jump);
   imove = find(xRef < 0);
   xRef(imove) = xRef(imove) + 2*pi;
   imove = find(xRef >= 2*pi);
   xRef(imove) = xRef(imove) - 2*pi;
   clear imove pp qq indx
end
%
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%==== Construct Sine coeff. for function expansion (size M-1). ========
%     sF(n) is n_th Sine Fourier Coeff. of fct.
%
Temp = fft(fct);
sF   = -imag(Temp([2:Np+1]))/M;
flag = max(abs(sF));
if flag <= 1e-7
   sF = 0*sF;
end
%
%==== Construct Cosine coeff. for function expansion (size M-1). =============
%     cF(n) is n_th Cosine Fourier Coeff. of fct.
%
cF   = real(Temp([2:Np+1]))/M;
cF_0 = real(Temp(1))/N;
flag = max(abs([cF, cF_0]));
if flag <= 1e-7
   cF = 0*cF; cF_0 = 0;
end
%
%==== Plot exact function. ===================================================
%
Maxf = max(fct);
Minf = min(fct);
range = Maxf - Minf + 1e-10;
Maxf = Maxf + 0.2*range;
Minf = Minf - 0.2*range;
Nrange = num2str(range);
fprintf(['\n Errors are relative to [max(F(x)) - min(F(x))] = ', ...
        Nrange, ' \n'])
figure
whitebg('k')
plot(x, fct, '-m', LW, 3)
axis([0 2*pi Minf Maxf])
ylabel('F(x)', FN, HNB, FS, 18)
xlabel('x', FN, HNB, FS, 18)
title('Exact function.', FN, HNB, FS, 20)
grid on
fprintf('\n Press any key to continue.\n')
pause
%==== Plot Fourier Coefficients. =============================================
figure
whitebg('k')
plot(sF, 'oc', MS, 6)
grid on
hold on
plot(sF, '-c')
ylabel('sF_n', FN, HNB, FS, 18)
xlabel('Fourier index n.', FN, HNB, FS, 18)
title('Sine Fourier Coefficients.', FN, HNB, FS, 20)
pause
nnn=0:1:Np;
%
figure
whitebg('k')
plot(nnn, [cF_0, cF], 'oc', MS, 6)
grid on
hold on
plot(nnn, [cF_0, cF], '-c')
ylabel('cF_n', FN, HNB, FS, 18)
xlabel('Fourier index n.', FN, HNB, FS, 18)
title('Cosine Fourier Coefficients.', FN, HNB, FS, 20)
pause
%
figure
whitebg('k')
power = sqrt([cF_0, cF].^2 + [0, sF].^2);
power = power + 1e-20;  % This is so 0's points in power are not dropped
                        % from the log plots. Error is bigger than this,
                        % anyway!
semilogy(nnn, power, '-r', LW, 3)
grid on
title('Semilog plot of Power Spectrum.',  FN, HNB, FS, 20)
ylabel('Power Spectrum.', FN, HNB, FS, 18)
xlabel('Fourier index n.',  FN, HNB, FS, 18)
pause
%
figure
whitebg('k')
loglog(nnn, power, '-r', LW, 3)
grid on
title('Loglog plot of Power Spectrum.',  FN, HNB, FS, 20)
ylabel('Power Spectrum.', FN, HNB, FS, 18)
xlabel('Fourier index n.',  FN, HNB, FS, 18)
pause
%
%==== Reconstruction and plot. ===============================================
%
rec = cF_0 * ones(1,N);
if Ndis >= 1
   recRef = cF_0 * ones(size(xRef));
end
Temp = abs(rec - fct);
goof = (max(Temp))/range;
%
figure
whitebg('k')
plot(x, rec,  '-r', LW, 3)
axis([0 2*pi Minf Maxf])
hold on
grid on
plot(x, fct)
xlabel('x', FN, HNB, FS, 18)
ylabel('y = F', FN, HNB, FS, 18)
title('Function and Mean Value.', FN, HNB, FS, 20)
pause
hold off
n = 1;
while n < Np+1
   Ndiff = min(Np-n+1, Nskip);
   nT=n:(Ndiff+n-1); nT = nT';
   recTemp    = diag(sF(nT))*sin(nT*x)    + diag(cF(nT))*cos(nT*x);
   if Ndis >= 1
      recTempRef = diag(sF(nT))*sin(nT*xRef) + diag(cF(nT))*cos(nT*xRef);
   end
   %  This here is a little trick to avoid putting sine and cosine evaluations
   %  inside a for loop.   Reason is sine and cosines are very slow in MatLab,
   %  so it is best to "vectorize" their evaluation to get it all done at once
   %  --- rather than inside a loop.  Else the first statement inside the loop
   %  should be replaced by:
   %                       rec = rec + sF(n)*sin(n*x) + cF(n)*cos(n*x); n=n+1;
   %  A real killer if Nskip and N are large!!
   for nn=1:Ndiff
       rec    = rec    + recTemp(nn, :);
       if Ndis >= 1  
          recRef = recRef + recTempRef(nn, :);
       end 
       n=n+1;
       Temp = abs(rec - fct);           % DO NOT CALCULATE ERROR THIS WAY AT
       Temp(Jump) = zeros(size(Jump));  % DISCONTINUITIES, because there are
       goofn = (max(Temp))/range;       % oscillations. Jump skips them.
       %
       if Ndis >= 1
          for jj=1:Ndis
              Maxdisrec(jj) = max(recRef((nwidth*(jj-1)+1):nwidth*jj));
              Mindisrec(jj) = min(recRef((nwidth*(jj-1)+1):nwidth*jj));
          end
          goofn = max(max(abs(Maxdisrec - Maxdisf))/range, goofn);
          goofn = max(max(abs(Mindisrec - Mindisf))/range, goofn);
       end
       %
       goof = [goof, goofn];
   end
   goofT = num2str(goofn);
   plot(x, rec,  '-r', LW, 3)
   axis([0 2*pi Minf Maxf])
   hold on
   grid on
   plot(x, fct)
   TEXT = num2str(n-1);
   title(['Fourier series added up to N = ', TEXT, '.'], FN, HNB, FS, 20)
   xlabel(['Relative error in approximation = ', goofT], FN, HNB, FS, 18)
   ylabel('y = F', FN, HNB, FS, 18)
   pause
   hold off
end
%
%==== Plot error now. =================================================
%
[n nn] = size(goof);
nnn = 0:nn-1;
figure
whitebg('k')
plot(nnn, goof,  '-r', 'LineWidth', 3)
grid on
title('Plot of relative error in approximation.',  FN, HNB, FS, 20)
ylabel('Error.', FN, HNB, FS, 18)
xlabel('Number of terms in Fourier sum.',  FN, HNB, FS, 18)
pause
%
figure
whitebg('k')
semilogy(nnn, goof,  '-r', LW, 3)
grid on
title('Semilog plot of relative error in approximation.',  FN, HNB, FS, 20)
ylabel('Error.', FN, HNB, FS, 18)
xlabel('Number of terms in Fourier sum.',  FN, HNB, FS, 18)
pause
%
figure
whitebg('k')
loglog(nnn, goof,  '-r', LW, 3)
grid on
title('Loglog plot of relative error in approximation.',  FN, HNB, FS, 20)
ylabel('Error.', FN, HNB, FS, 18)
xlabel('Number of terms in Fourier sum.',  FN, HNB, FS, 18)
pause
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n To clear all figures and continue, enter FLAG = 1.')
fprintf('\n To STOP NOW ....................., enter FLAG = 0.\n')
myFLAG = input('                                          FLAG = '); 
if myFLAG > 0.5
   close all
end
end
clear  Maxf      TEXT      nT        Temp      flag      nn        CLR
clear  yb        FN        Minf      VERSION   nnn       yt        FS
clear  N         goofT     HNB       NN        h         range     jj
clear  LW        h_fig     M         Ndiff     OPTION    dyt       goofn
clear  jn        recTemp   MM        OPT       dyb       myFLAG    MS
clear  n         smooth    DOTS      TEXTT     YT        dl        dy
clear  nop       x_m       y_m       yl        BGC       Nrange    recTempRef
%
% CLEAN UP xRef so it is ordered:
%
if Ndis >= 1
   [xRef, II] = sort(xRef);
   recRef = recRef(II);
end
clear II Ndis Maxdisf Mindisf Maxdisrec Mindisrec
%
% Ndis       number of discontinuities.
% Maxdisf    maximum values at discontinuities of function.
% Mindisf    minimum values at discontinuities of function.
% Maxdisrec  maximum values near discontinuities of approximation.
% Maxdisrec  minimum values near discontinuities of approximation.
%
% EOF