% Using SVD to gain insight into an inverse problem
% here: deblurring a 1-d signal

n = 100;
w = 3;

v1 = ones(n,1);  % vector of all ones

X = zeros(n,n);       % X is the  blur operator, 
                      % here it's (2w+1)-point moving average with wraparound
for i=-w:w,
   X = X + spdiags(v1,i,n,n);    % adding near-diagonal elements
end
   
corner = ones(w,w);
X(1:w,n-w+1:n) = triu(corner);
X(n-w+1:n,1:w) = tril(corner);   % adding 

X = X/(2*w+1);                % rescaling for row totals = 1; not strictly necessary
figure(1)
imagesc(X);   colormap gray

%============================ Blurring

m = [1:22 22*cos(0.1*(1:28)) zeros(1,15) (-10)*ones(1,10) 10*sin((1:25)/2)  ]';  % some made-up signal
plot(m)

d = X*m;
   hold on
   plot(d,'*r')
   hold off
   
 [U,S,V] = svd(X,0);     % compact SVD
 figure(2)
  plot(diag(S))
 figure(3); bookfonts
  imagesc(U)      
    colormap gray
    title('Visualizing matrix U')
  
 
 figure(4) ; bookfonts
 plot(U(:,1)) 
   hold on 
   plot(U(:,10),'*b')
   plot(U(:,50),'.-r')
   title('First, 10th and 50th basis functions')
    hold off
    
 % ======================= Deblurring 

  d1 = d + randn(size(d));    % norm(d)=100  so this is about 1% noise 

  sigmas = diag(S);   
  m1 = V*((U'*d1)./sigmas);   % possibly a faster implementation of V * S^(-1) * U' * d
  
  figure(5); bookfonts
  plot(m1)
     title('Deblurring a noisy signal')
     
 min(sigmas)    % this provides a clue of what went wrong 
 
  k = 20;                     %  truncation level
   part1 = (U'*d1)./sigmas;
   part1(k+1:end) = 0;         % simply truncate all the "higher order" information
   m2 = V*part1;
   
   figure(6); bookfonts
   plot(m)
     hold on 
     plot(m2,'*r')
      title(['Original and restored signal, using k = ' num2str(k) ' basis vectors']); 
     hold off
     
     % the conclusion is: some degree of smoothing is inevitable, we need to
     % pick k as large as possible while making sure we do not "amplify noise"
     
     % -- this is sometimes known as "bias/variance tradeoff"