%% Matlab or Octave script
%% Who:Abdallah Sayyed-Ahmad
%% What: Least Square approximation (polynomials)

% get some data points from a known function, but also can be obtained from measurment data
x=[-4:0.5:4];
f=cos(x);
N=numel(x);
M=16;%order of the polynoial approximation
X=zeros(N,M+1);

%compute the Vandermond martix
for  i=1:N
  for j=1:M+1
    X(i,j)=x(i)^(j-1);
    endfor
endfor

%solve X'Xa=X'b
a=X\f';
% evaluate the solution for the range of data at a higher resolution
xx=[-4:0.05:4];
NN=numel(xx);
yy=zeros(1,NN);
for  i=1:NN
  yy(i)=0;
  for j=1:M+1
    yy(i)=yy(i)+a(j)*xx(i)^(j-1);
    endfor
endfor
%plot least square approximation solution
plot(xx,yy)
hold on
%plot original datapoints
scatter(x,f)