% 1D interpolation implemented with least-squares

clear all;

% 100 m's that we want to know
Nx=101;
dx=1.0;
xmin=0.0;
xmax = xmin+dx*(Nx-1);
x = xmin+dx*[0:Nx-1]';

% but we only have measured 40 of them (randomly chosen)
Nd=40;
rowindex=sort(unidrnd(Nx,Nd,1));
xd=xmin+dx*rowindex;

% rhe true data are a sine curve
f=0.02;
y = sin(2*pi*f*xd);

% plot the observed data
figure(1);
clf;
plot(xd,y,'r+');
hold on;
axis([xmin xmax -2 2]);
xlabel('x');
ylabel('y');

M=Nx;
m=zeros(M,1);

% top part of matrix is Gm=d
Nc=Nx;
N=Nd+Nc;
G=zeros(N,M);
d=zeros(N,1);
for p = [1:Nd]
    G(p,rowindex(p)) = 1.0;
    d(p)=y(p);
end

% bottom part pf matric is epsilon Dm = 0
epsilon = 1e-6;
a=epsilon*(1.0);
b=epsilon*(-2.0);
c=epsilon*(1.0);
f=epsilon*(-1.0);

% second derivative in interior
for p = [1:Nc-2]
    G(p+Nd,p) = a;
    G(p+Nd,p+1) = b;
    G(p+Nd,p+2) = c;
    d(p+Nd)=0.0;
end

% first derivative at left end
G(N-1,1)=a;
G(N-1,2)=f;
d(N-1)=0;

% first derivative at right end
G(N,M-1)=a;
G(N,M)=f;
d(N)=0;

% solve by least squared
m= inv(G'*G)*G'*d;
plot(x,m,'b');
hold on;

% evaluate error
e = y-m(rowindex(:));
E=e'*e;

title( sprintf('error = %f', E ));

