% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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