


clear
close all

% Illustration of some of the fundamental concepts of linear inverse
% problems with the example of fitting a straight line to some noisy data

% The data are supposed to mimick a reflection arrival 

% I. Generate some noise data along an x axis in the interval [0, 1]

% the x axis
n=5;
dx=2;
x=dx:dx:n*dx;

% True solution of the 2-parameter problem
% Intercept time 1 s, slope 1/3 s/km
m1t=1.;
m2t=1/3. ;

y=m1t+m2t*x;

% Plot true solution without noise

figure
plot(x,y,'x')
xlabel('x')
ylabel(' Data ')
title(' True noise-free data ')
hold on
plot(x,m1t+m2t.*x)
hold off

figure
% Generate noise 
pert=.5;
y=y+(rand(1,length(y))-.5)*2*pert;
plot(x,y,'x')
xlabel('x')
hold on
plot(x,m1t+m2t.*x)
ylabel(' Noisy Data ')
title(' True noisy data ')
hold off


%%%%%%%%%%%%%%%%%%% The linear inverse problem %%%%%%%%%%%%%%%%%%%%%%%%%

% the data vector
% we just rename to stick with the consistent formulation of G m = d
d = y;

% the "model" vector m (is two dimensions) 
m=zeros(1,2)';

% the system matrix G (the kernel, the forward operator, ...)
G=[ones(1,length(x));  x]';

% let us analyze G 
rank(G);

% Is it column rank (2)?

% Finally the solution to the overdetermined linear inverse problem

% for problem
d=[  1.4102    2.7626    2.8500    3.3633    4.0844];
m=inv(G'*G)*G'*d';


%%%%% plot the solution
figure
plot(x,d,'x')
xlabel('x')
hold on
plot(x,m(1)+m(2).*x)
ylabel(' Noisy Data ')
title(' Data and solution ')
hold off

disp(sprintf(' Inverted Slope: %g km/s, true: %g km/s',1/m(2),1/m2t));
disp(sprintf(' Inverted t_i: %g s, true: %g s',m(1),m1t)) 