clear all
format long

g = @(x) exp(-x).*cos(4*x);     % define the function to approximate

n = 3;                          % use 2n+1 exponentials
n_data_pts = 20;                % number of data points to approximate

fprintf('using %d terms in the approximation\n',2*n+1);

x = linspace(-1,1,n_data_pts)'; % the x-locations of the data points
y = g(x);                       % the y-locations of the data points


% Assemble the matrix A, a_ij = e^(k*x_i) for k=j-(n+1)
k_vector=-n:n;
for i=1:n_data_pts              % loop over the data points
    A(i,:) = exp(k_vector*x(i)); % assemble row i of the matrix
end

b = A\y                         % find the best fit
% b_pinv = pinv(A)*y            % get the same answer using pseudoinverse
residual_norm = norm(y-A*b);    % compute the norm of the residual
fprintf('||r(b)||_2/||b||_2 = %d\n\n',residual_norm/norm(b));

xx = linspace(-1,1);            % plot g and the best fit
yy = g(xx);                     % evaluate the true function at the plotting pts

% accumulate the best fit data with a loop
best_fit_yy = 0*xx;
for j=1:(2*n+1)
    k = j - (n+1);
    best_fit_yy = best_fit_yy + b(j)*exp(k*xx);
end

% plot g(x)
plot(xx,yy,'b','linewidth',2);
hold on
plot(x,y,'bx','markersize',12,'linewidth',2)
% plot the best fit
plot(xx,best_fit_yy,'k','linewidth',2);
legend('e^{-x}cos(4x)','samples','best fit')
title([num2str(2*n+1),' terms in approximation'],'fontsize',16)
hold off
