 close all; clear all;
%%-------------------- function to produce signal
 ff=inline('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 
%%--------------------smoothing/blurring matrix
%% 126x126 matrix 
 load('Blurmat');
%%--------------------Apply blurring to data
 zz = A*yy;
%%-------------------- add some noise to result 
 noiseLevel = 1.e-07;
 y1 = zz + noiseLevel*randn(m,1); 
%%-------------------- plot modified signal. Should
%%                     be close to original signal
% figure(2) 
 plot(xx,y1,'r--','linewidth',2) 
 title('Transformed data + noise','fontsize',18) 
 pause (2)
%%--------------------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 (2)
%%-------------------- recover by truncated SVD -- 
  yy2 = pinv(A,noiseLevel)*y1;
 %%yy2 = pinv(A,0.05)*y1;
 figure(4) 
 plot(xx,yy,'o',xx,yy2,'red') 
 title(' Recovery: TSVD','fontsize',18) 
 hold off;


