 close all; clear all;
%%-------------------- function to produce signal
 ff=@(x) sin(3*pi*x)-.2*cos(pi*5*x)+sin(pi*x)
 xx = [0:0.05:2*pi]'; 
 yy = ff(xx) ;
 m = length(xx); 
%%-------------------- plot orig. data
 figure(1)
 plot(xx,yy,'linewidth',2) 
 hold on
 title(' Original data','fontsize',18)   
 pause 
%%-------------------- load blurring matrix - 
%% 126x126 matrix 
 load('Blurmat');
%%-------------------- transform data by A 
 zz = A*yy;
%%-------------------- add some noise to result 
 noiseLevel = 0.000001;
 y1 = zz + noiseLevel*randn(m,1); 
%%-------------------- plot modified signal -- 
%% figure(2) 
 plot(xx,y1,'r--','linewidth',2) 
 title('Transformed data + noise','fontsize',18) 
 pause
%%--------------------recover blurred & perturbed
%%                    signal by EXACT INVERSION:
 yy1= A \ y1;   
 figure(3);
 plot(xx,yy,'o',xx,yy1,'red') 
 title('Recovery: direct inversion','fontsize',18)
 pause 
%%-------------------- recover by truncated SVD
   yy2 = pinv(A,noiseLevel)*y1;
%%  yy2 = (A+1.e-05*eye(size(A)))\y1;

 %%yy2 = pinv(A,0.05)*y1;
 figure(4) 
 plot(xx,yy,'o',xx,yy2,'red') 
 title(' Recovery: TSVD','fontsize',18) 
 hold off;


