% 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))
  
 d1 = d + randn(size(d));    % norm(m)=100 to this is about 1% noise 
 
 % ======================= Deblurring