%MATLAB script to generate a figure like Figure 6.7 on page 271 of conte/deboor
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
% generate smooth 2pi-periodic function roughly like the one in that figure
x = (-1:.2:1)*pi;
y = [0.7930   -0.2216   -0.5598   -0.3557    0.2741    0.6006    0.3149 ...
     0.3207    0.7405    1.0738    0.7930];
   % in order to get a 2*pi-periodic curve, extend the data beyond the interval
   % [-pi .. pi], spline interpolate, then restrict the resulting spline to
   % the interval [-pi .. pi]
X = [x(6:10)-2*pi, x, x(2:6)+2*pi]; Y = y([6:10 1:11 2:6]);

sp = CALCCF(X,spline_cd(X,[Y;zeros(1,length(Y))]));
    % look at the data to see that it makes sense to choose zero slope
    % at  X(1)  and  X(end) .

pp = struct('breaks',x,'coefs',sp.coefs(:,6:15));

% check periodicity
h = pp.breaks(end)-pp.breaks(end-1);
fprintf('the difference between the two endslopes is %12.5e\n', ...
pp.coefs(2,1) - (pp.coefs(2,end)+h*(2*pp.coefs(3,end)+h*3*pp.coefs(4,end))))
fprintf('the difference between the two second derivs is %12.5e\n', ...
2*pp.coefs(3,1) - (2*pp.coefs(3,end)+h*6*pp.coefs(4,end)))

%  generate a noisy version of this function
xx =(-60:59)/60*pi; yy = PCUBIC(pp,xx); 
%  make this reproducible by setting the random number generator to its
%  default value.
   rng('default'), yyn = yy+(rand(1,120)-.5)/20;

plot(xx,yyn,'.'), title('noisy data'), pause

%  compute the discrete Fourier transform of the noisy data ...
[Z1,Z2,out] = fft_cd( yyn,yyn,1);
%  ... and plot the absolute values of the resulting discrete Fourier
%  coefficients
eval(['zz = Z',num2str(out),';'])
xxx = repmat(1:120,3,1); azz = abs(zz(:)'); 
zzz = [zeros(size(azz));azz;repmat(NaN,size(azz))];
figure, plot(xxx(:),zzz(:),'-')
title('absolute values of discrete Fourier coefficients'), pause

%  evaluate the corresponding Fourier series at the  xx ...
[Z3,Z4,outs] = fft_cd(conj(zz/120),yyn,1);
%  ... and compare with the original data; in exact arithmetic, their
%  difference should be zero.
eval(['zzr = Z',num2str(2+outs),';'])
fprintf('error in recovered values is %10.3e\n',max(abs(yyn(:)-zzr(:))))

%  obtain smoothed values at the  xx  by setting to zero all Fourier
%  coefficients smaller in absolute value than 1 percent of the absolutely
%  largest Fourier coefficient.

zz(azz<max(azz(:))/100) = 0;
[Z3,Z4,outs] = fft_cd(conj(zz/120),yyn,1);
eval(['zzs = Z',num2str(2+outs),';'])
figure, plot(xx,real(zzs(:)'),'r-'), hold on, plot(xx,yyn,'k.')
title('smoothed values from noisy data')
