% Do regularised inversion using Truncated SVD
%
n = 64;
sigma = 0.1; % std.dev. of added white noise
f = zeros(n,1);
f(12) = -1; f(20:28) = 1.5; f(32:36) = 2; % arbitrary function
[y,K] = linblur(f,0.04);
%[y,K] = weighted_linblur(f,0.04);
yn = y + sigma*randn(n,1); % add white noise, 

figure (1);clf; hold on;
plot(y);
plot(yn,'g');
plot(f,'r');

[U,W,V] = svd(K);
disp(['det K : ',num2str(det(K))]);
disp(['condition number : ',num2str(W(1,1)/W(n,n))]);
figure(2); 
subplot(2,3,1);plot(diag(W));
subplot(2,3,2);plot(U(:,1));
subplot(2,3,3);plot(U(:,2));
subplot(2,3,4);plot(U(:,4));
subplot(2,3,5);plot(U(:,8));
subplot(2,3,6);plot(U(:,16));

figure(3);clf; hold on;
plot(f,'r');
figure(1);clf; hold on;
plot(y,'b'); plot(yn,'g');
figure(5); clf; imagesc(K);
%%
fi = zeros(size(f));
for k=1:floor(0.45*n)
    fi = fi + U(:,k)'*yn *V(:,k)/W(k,k);
    figure(3);
    plot(fi,'k');
    pause(1);
    figure(1);
    plot(K*fi,'k');
    pause(2);
end
figure(4);
fp = K\yn;
plot(fp,'m');


%% produce a "publication" figure;


fi4 = zeros(size(f)); fi16 = fi4; fi32 = fi4;
for k = 1:4
fi4 = fi4 + U(:,k)'*yn *V(:,k)/W(k,k);
end
yi4 = K*fi4;
for k = 1:16
fi16 = fi16 + U(:,k)'*yn *V(:,k)/W(k,k);
end
yi16 = K*fi16;

for k = 1:32
fi32 = fi32 + U(:,k)'*yn *V(:,k)/W(k,k);
end
yi32 = K*fi32;

figure(10); clf;
subplot(1,2,1);
plot(y); hold on;
plot(yn,'g');
plot(yi4,'k');
plot(yi16,'-k+');
plot(yi32,'-.k.');

subplot(1,2,2);
plot(f,'r'); hold on;
plot(fi4,'k');
plot(fi16,'-k+');
plot(fi32,'-.k');

%% and with Tikhonov Filtering

for alpha = W(1,1):-W(1,1)/20:W(1,1)/40
    disp(['\alpha = ',num2str(alpha)]);
    fi = zeros(size(f));
for k=1:n
    fi = fi + U(:,k)'*yn *V(:,k)*W(k,k)/(W(k,k)^2 + alpha);
end
    figure(3);
    plot(fi,'m');
    pause(1);
    figure(1);
    plot(K*fi,'m');
    pause(2);
end
figure(4);hold on;
fp = K\yn;
plot(fi,'--c');


%% produce a "publication" figure for TK0;

alpha1 = 0.9*W(1,1);
alpha2 = 0.2*W(1,1);
alpha3 = 0.04*W(1,1);
fa1 = zeros(size(f)); fa2 = fa1; fa3 = fa1;
for k=1:n
    fa1 = fa1 + U(:,k)'*yn *V(:,k)*W(k,k)/(W(k,k)^2 + alpha1);
    fa2 = fa2 + U(:,k)'*yn *V(:,k)*W(k,k)/(W(k,k)^2 + alpha2);
    fa3 = fa3 + U(:,k)'*yn *V(:,k)*W(k,k)/(W(k,k)^2 + alpha3);
end


ya1 = K*fa1;
ya2 = K*fa2;
ya3 = K*fa3;


figure(11); clf;
subplot(1,2,1);
plot(y); hold on;
plot(yn,'g');
plot(ya1,'k');
plot(ya2,'-k+');
plot(ya3,'-.k.');

subplot(1,2,2);
plot(f,'r'); hold on;
plot(fa1,'k');
plot(fa2,'-k+');
plot(fa3,'-.k');