clear

n = 11
x0 = ones(n,1)
x1 = linspace(-1,1,11)'
x2 = x1.*x1
figure(1)
clf
plot(x1,[x0 x1 x2])
hold on
plot(x1,[x0 x1 x2], '*')
grid on
xlabel('x')
ylabel('f(x)')
title('not normalized')

% norm is just the length of a vector - in 2D l = sqrt( x^2 + y^2 )
norm(x0)
norm(x1)
norm(x2)

dot(x0,x1)
dot(x0,x2)
dot(x1,x2)

x0n = x0/norm(x0)
x1n = x1/norm(x1)
x2n = x2/norm(x2)

figure(2)
clf
plot(x1,[x0n x1n x2n])
hold on
plot(x1,[x0n x1n x2n], '*')
grid on
xlabel('x')
ylabel('f(x)')
title('normalized')

dot(x0n,x1n)
dot(x0n,x2n)
mean(x2n)
x2n=x2n-mean(x2n)

dot(x0n,x2n)
dot(x1n,x2n)
dot(x0n,x1n)
% Another way of writing the dot product or scalar product 
% or inner product
x0n'*x1n

% replot with orthogonal 2nd power
figure(2)
clf
plot(x1,[x0n x1n x2n])
hold on
plot(x1,[x0n x1n x2n], '*')
grid on
xlabel('x')
ylabel('f(x)')
title('normalized')

% check orthonality of a whole set of vectors at once

H = [x2n x1n x0n]

H'*H

% orthonormalize
x2n = x2n/norm(x2n)
H = [x2n x1n x0n]
H'*H

% but not row orthogonal
H*H'
%[u,v,w] = svd(H)
%H - u*v*w'

%a = u'*u-eye(size(u,1));
%max(abs(a(:)))

%a = w'*w-eye(size(w,1));
%max(abs(a(:)))

% make a simulate, noise-free measurement
g = 4*x2n - 2*x1n + 1*x0n

% can get values back by projections onto a vector
dot(g,x2n)
dot(g,x1n)
dot(g,x0n)

% Let's add one more power
g2 = g+rand(size(g))
g2-g
x3n=x1n.^3
norm(x3n)
x3n=x3n/norm(x3n)
H2 = [x3n x2n x1n x0n]
H2'*H2

plot(x1,x3n,'g')
plot(x1,x3n,'g^')

% but not orthogonal
g2 = g - 6*x3n
dot(g2,x0n)
dot(g2,x1n)
dot(g2,x2n)
dot(g2,x3n)

% the fit coefficients change with the order of the fit
f = H2\g2
f = H\g2

%g2n = g2+(rand(size(g2))-.5)*.05
%f = H2\g2n
%f = H\g2n

% Gram-Schmidt orthogonalization
x3n = x3n - x1n*dot(x1n,x3n)/norm(x1n);
x3n = x3n/norm(x3n)
H2(:,1) = x3n
H2'*H2

% replot all 4 vectors with orthogonal 2nd power
figure(2)
clf
plot(x1,[x0n x1n x2n x3n])
hold on
plot(x1,[x0n x1n x2n, x3n], '*')
grid on
xlabel('x')
ylabel('f(x)')
title('normalized')

% Make a new version of sampled data that uses the orthogonal
% 3rd power
g2 = g - 6*x3n
dot(g2,x0n)
dot(g2,x1n)
dot(g2,x2n)
dot(g2,x3n)

% now fit does not change - that's why an orthogonal or even better,
% orthonormal basis is useful
f = H2\g2
f = H\g2