function decmalp_a();
%            Deconvolution is an ill posed inverse problem
% f(t) input signal
% g(t) output signal 
% g(t)=h(t)*f(t)+e(t)
% We also have the respose of the system u(t) to a step function
% u(t)=\int h(t) dt   step function response
%
% Author: Ali Mohammad-Djafari
% Revision date: 10/1/2006
%
% 1. Simulate the data
% 2. No noise, everything works
% 3. With noise, no simple method work
% 4. Ad hoc methods
% 5. Wienner Filtering
% 6. Regularization (zero order)
% 7. Regularization (first order)
% 8. First estimate h from u, then f using g and h
% 9. Blind deconvolution: estimate h and f from y
%10. Blind deconvolution: Implementation by FFT
% 0. Quit

clear all, clc, help decmalp_a
figure(1);h=gcf;set(h,'Position',[400,400,500,400]);

while 1
 %clc;help decmalp_a
 %nm1 = input('Votre choix? =>');
 %if (nm1 <= 0), break, end

 nm1=menu('Deconvolution is an ill posed inverse problem','Simulate the data', 'No noise: all methods work','With noise: no simple methods work','Ad hoc methods','Wienner filtering','Regularization (zero order)','Regularization (first order)', 'First estimate h from u and then f from g and h','Blind Deconvolutio: estimate h and f from g','Blind Deconvolution: Implementation using FFT','Quit');
%---------------------------------------------------------------------------
 if (nm1 == 1)
  disp('Generate data')
  nm2=menu('Choose one of these examples','Example1','Example 2','Example 3','Example 4');
  %nm2=1;
  switch nm2
  %-------------------------
  case 1
   t=[0:31]';
   f=[t;zeros(31,1)];f=(1/max(f))*f;lf=length(f);
   lh=32;a=[1;-1.6;.7];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);
   f=conv(f,h);f=f(3:lf+2);f=f/max(f);
 
   lh=32;a=[1;-1.2;.3];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);

   g=conv(f,h);g=g(1:lf);lg=length(g);
   h1=[h;zeros(lf-lh,1)];lh1=length(h1);
   u=cumsum(h);lu=length(u);
   u1=cumsum(h1);lu1=length(u1);
  %-------------------------
  case 2
   t=[0:31]';
   f=[t;zeros(31,1)];f=(1/max(f))*f;lf=length(f);
   lh=32;a=[1;-1.6;.7];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);
   f=conv(f,h);f=f(3:lf+2);f=f/max(f);

   t=[0:31]';h=.1*exp(-.1*t).*cos(.1*t);lh=length(h);
   %lh=32;a=[1;-1.1;.2];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);

   g=conv(f,h);g=g(1:lf);lg=length(g);
   h1=[h;zeros(lf-lh,1)];lh1=length(h1);
   u=cumsum(h);lu=length(u);
   u1=cumsum(h1);lu1=length(u1);
   %-------------------------
   case 3
   f=[zeros(9,1);1;zeros(10,1);2;zeros(11,1);];
   f=.5*[f;zeros(4,1);1;zeros(13,1);1.5;zeros(12,1)];lf=length(f);

   t=[0:31]';h=.1*exp(-.1*t).*cos(.1*t);lh=length(h);

   g=conv(f,h);g=g(1:lf);lg=length(g);
   h1=[h;zeros(lf-lh,1)];lh1=length(h1);
   u=cumsum(h);lu=length(u);
   u1=cumsum(h1);lu1=length(u1);
   %-------------------------
   case 4
   f=[zeros(5,1);ones(5,1);zeros(4,1);2*ones(4,1);zeros(6,1);1*ones(6,1);];
   f=[f;zeros(4,1);3*ones(4,1);zeros(6,1);1*ones(6,1)];f=[f;f];lf=length(f);
 
   t=[0:31]';h=.1*exp(-.1*t).*cos(.1*t);lh=length(h);

   g=conv(f,h);g=g(1:lf);lg=length(g);
   h1=[h;zeros(lf-lh,1)];lh1=length(h1);
   u=cumsum(h);lu=length(u);
   u1=cumsum(h1);lu1=length(u1);
  end
  disp('Save data in fuhg.mat')
  save fuhg.mat f u h g h1 u1

  figure(1),clf;plot(f,'r'),hold on;plot(u1,'b'),plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  th=text(33,.9,'f(t)');set(th,'Color','r')
  th=text(20,1.1,'u(t)');set(th,'Color','b')
  th=text(5,.9,'h(t)');set(th,'Color','r')
  th=text(40,.7,'g(t)');set(th,'Color','b')
  hold off
  print -depsc Eps/fuhg.eps;
  %disp('Appughz sur une touchh pour continuer'), pause

  disp('Add noise')
  gb=g+.01*randn(length(g),1);
  ub=u+.01*randn(length(u),1);
  u1b=u1+.01*randn(length(u1),1);
  disp('Save noisy data in fuhgb.mat')
  save fuhgb.mat  u h g h1 u1 gb ub u1b 

  figure(1),clf;plot(f,'r'),hold on;plot(u1,'b'),plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  th=text(33,.9,'f(t)');set(th,'Color','r')
  th=text(20,1.1,'u(t)');set(th,'Color','b')
  th=text(5,.9,'h(t)');set(th,'Color','r')
  th=text(40,.7,'g(t)');set(th,'Color','b')
  hold off,pause(.1)
  print -depsc Eps/fuhgb.eps;
 end
%---------------------------------------------------------------------------
 if (nm1 == 2)
  disp('No noise : Any method may work')
  disp('Load data without noise')
  load fuhg.mat
  disp('Without noise even simple methods works. For example:')
  disp('  -- Just compute the derivative of u(t) to obtain h(t)')
  disp('  -- Use Inverse Filtering to do deconvolution.')

  lf=length(f);lg=length(g);lh=length(h);lu=length(u);
  lh1=length(h1);lu1=length(u1);

  hh=[0;u1(2:lu1)-u1(1:lu1-1)];
  H=fft(hh); G=fft(g); Fh=G./(H+.1); 

  fh=real(ifft(Fh));

  %Gh=H.*fft(fh);gh=real(ifft(Gh));
  gh=conv(fh,hh);gh=gh(1:lf);

  uh=cumsum(hh)+u(1);

  figure(1),clf;plot(f,'r'),hold on;plot(u1,'b'),plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  
  plot(uh,'+g');plot(uh,'g');th1=text(20,1.1,'u(t)');set(th1,'Color','b');th2=text(25,1.1,'uh(t)');set(th2,'Color','g')
  plot(8*hh,'+g');plot(8*hh,'c');th1=text(5,.9,'h(t)');set(th1,'Color','r'),th2=text(10,.9,'hh(t)');set(th2,'Color','g')
  plot(fh,'+c');plot(fh,'c');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','c')
  plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
  hold off 
  eval(['print -depsc Eps/xuhgh',num2str(nm1),'.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 > 2)
  disp('Load the noisy data')
  load fuhgb.mat
  g=gb;
  u=ub;u1=u1b;
  lf=length(f);lg=length(g);lh=length(h);lu=length(u);
  lh1=length(h1);lu1=length(u1);
 end

 if (nm1 >= 3)&(nm1 <= 6)
  hh=[0;u1(2:lu1)-u1(1:lu1-1)];H=fft(hh); 
 end
%---------------------------------------------------------------------------
 if (nm1 == 3)
  disp('With noise, no naive methods work')
  G=fft(g); Fh=G./H; fh=real(ifft(Fh));
  gh=conv(fh,hh);gh=gh(1:lf);
  uh=cumsum(hh)+u(1);
  
  figure(1),clf;plot(f,'r'),hold on;plot(u1,'b'),plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  plot(uh,'+g');plot(uh,'g');th1=text(20,1.1,'u(t)');set(th1,'Color','b');th2=text(25,1.1,'uh(t)');set(th2,'Color','g')
  plot(8*hh,'+g');plot(8*hh,'c');th1=text(5,.9,'h(t)');set(th1,'Color','r'),th2=text(10,.9,'hh(t)');set(th2,'Color','g')
  plot(fh,'+c');plot(fh,'c');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','c')
  plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
  hold off
  eval(['print -depsc Eps/xuhgh',num2str(nm1),'.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 == 4)
  disp('Ad hoc methods to avoid DIVIDE BY ZERO') 
  eps1=.1;eps1=getnew(eps1,'Seuil');
  iz=find(abs(H)<eps1);H(iz)=eps1*ones(size(iz));
  G=fft(g); Fh=G./H;
  fh=real(ifft(Fh));
  gh=conv(fh,hh);gh=gh(1:lf);

  figure(1),clf;plot(f,'r'),hold on;plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  plot(fh,'+g');plot(fh,'g');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','g')
  plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
  hold off
  eval(['print -depsc Eps/xuhgh',num2str(nm1),num2str(10*eps1),'.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 == 5)
  disp('Wienner Filtering') 
  r=20;r=getnew(r,'SNR in dB');  %r=10 20 30
  r1=10^(-r/20);
  G=fft(g); 
  Fh=G.*conj(H)./(abs(H).^2+r1); 
  fh=real(ifft(Fh));
  gh=conv(fh,hh);gh=gh(1:lf);

  figure(1),clf;plot(f,'r'),hold on;plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  plot(fh,'+g');plot(fh,'g');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','g')
  plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
  hold off
  eval(['print -depsc Eps/xuhgh',num2str(nm1),num2str(r/10),'.eps']);
 end
%---------------------------------------------------------------------------
 if ((nm1 >= 6) & (nm1 < 8))
  hh1=[0;u1(2:lu1)-u1(1:lu1-1)]; 
  H=toeplitz([zeros(lg,1)],hh1);H=H'; 
  HtH=H'*H;
  if (nm1 == 6)
   disp('Regularization (zero order): Try the values of lambda .1 1 10')
   lambda=.1;lambda=getnew(lambda,'Regularization coefficient lambda');
   DtD=eye(lg,lg);
   fh=(HtH+lambda*DtD)\(H'*g);
  end
%---------------------------------------------------------------------------
  if (nm1 == 7)
   disp('Regularization (First order): Try the values of lambda .1 1 10')
   lambda=.1;lambda=getnew(lambda,'Regularization coefficient lambda');
   D=convmtx([1,-1],lg-1);DtD=D'*D;
   fh=(HtH+lambda*DtD)\(H'*g);
  end
  gh=conv(fh,hh);gh=gh(1:lf);

  figure(1),clf;plot(f,'r'),hold on;plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  th=text(40,.7,'g(t)');
  plot(fh,'+g');plot(fh,'g');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','g')
  plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
  hold off
  eval(['print -depsc Eps/xuhgh',num2str(nm1),num2str(100*lambda),'.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 == 8)
  disp('First estimate h(t) from u(t)')
  U=toeplitz([1;zeros(lu1-1,1)],ones(lu1,1));U=U'; 
  DtD=eye(lg,lg);
  lambda=1;lambda=getnew(lambda,'Regularization Coefficient for estimating h');
  hh=(U'*U+lambda*DtD)\(U'*u1);
   
  disp('Then estimate x(t) using estimated h')
  H=toeplitz([hh(1);zeros(lg-1,1)],hh);H=H';HtH=H'*H; 
  lambda=.02;lambda=getnew(lambda,'Regularization Coefficient for estimating x');
  DtD=eye(lg,lg);
  fh=(HtH+lambda*DtD)\(H'*g);
  gh=conv(fh,hh);gh=gh(1:lf);

  figure(1),clf;plot(f,'r'),hold on;plot(u1,'b'),plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
  th=text(20,1.1,'u(t)');set(th,'Color','b');
  plot(fh,'+g');plot(fh,'g');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','g')
  plot(8*hh,'+c');plot(8*hh,'c');th1=text(5,.9,'h(t)');set(th1,'Color','r'),th2=text(10,.9,'hh(t)');set(th2,'Color','c')
  plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
  hold off
  eval(['print -depsc Eps/xuhgh',num2str(nm1),'.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 == 9)
  disp('Blind Deconvolution: Estime h and f from g')
  lambda1=300;lambda1=getnew(lambda1,'Regularization Coeff for h');
  lambda2=1;lambda2=getnew(lambda2,'Regularization Coeff for f');
  fh=g;lh=31;J1=0;J2=0;dJ=-1;
 
  k=0;
  while((dJ<=0)&(abs(dJ)>.001))
   k=k+1;
   if(lambda1>.1),lambda1=.7*lambda1;end
   if(lambda2>.01),lambda2=.7*lambda2;end
   %A=toeplitz([fh(1);zeros(lf-1,1)],fh);A=A'; 
   A=toeplitz([fh(1);zeros(lh-1,1)],fh);A=A'; 
   D=convmtx([1,-1],lh-1);DtD=D'*D;
   %D=egh(lh,lh);DtD=D;
   hh=(A'*A+lambda1*DtD)\(A'*g); 
   hh=(1/sum(hh))*hh;
   iz=find(hh < 0);if(~isempty(iz)),hh(iz)=zeros(length(iz),1);end;
   hh1=[hh;zeros(lg-lh,1)];
   dhh=D*hh;
  
   %hh(lh+1:lf)=zeros(lf-lh,1);
  
   H=toeplitz([hh(1);zeros(lg-1,1)],hh1);H=H'; 
   D=convmtx([1,-1],lg-1);DtD=D'*D;
   %DtD=eye(lg,lg);
   fh=(H'*H+lambda2*DtD)\(H'*g);
   dfh=D*fh;
   gh=H*fh;
  
   dg=g-gh;
   J2=dg'*dg+lambda1*(dhh'*dhh)+lambda2*(dfh'*dfh);disp([k,J2])
  
   %J2=dg'*dg
   if(k>1),dJ=J2-J1;J1=J2;else,dJ=-1;J1=J2;end

   figure(1),clf;plot(f,'r'),hold on;plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
   plot(fh,'+g');plot(fh,'g');th1=text(32,.9,'f(t)');set(th1,'Color','r'),th2=text(38,.9,'fh(t)');set(th2,'Color','g')
   plot(8*hh,'+c');plot(8*hh,'c');th1=text(5,.9,'h(t)');set(th1,'Color','r'),th2=text(10,.9,'hh(t)');set(th2,'Color','c')
   plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
   hold off
   if(k==1),eval(['print -depsc Eps/xuhgh',num2str(nm1),'1.eps']);end
   pause(1)
  end
  eval(['print -depsc Eps/xuhgh',num2str(nm1),'f.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 == 10)
  disp('Blind Deconvolution: Implementation by FFT')
  lambda1=20;lambda1=getnew(lambda1,'Regularization Coeff for h');
  lambda2=.1;lambda2=getnew(lambda2,'Regularization Coeff for f');
  fh=g;lh=31;J1=0;J2=0;dJ=-1;
 
  Yw=fft(g);
  Dhw=fft([-1;1;zeros(lg-2,1)]);Dhw2=abs(Dhw).^2;
  Dfw=fft([-1;1;zeros(lg-2,1)]);Dfw2=abs(Dfw).^2;
  k=0;
  while((dJ<=0)&(abs(dJ)>.001))
   k=k+1;
   if(lambda1>.1),lambda1=.7*lambda1;end
   if(lambda2>.01),lambda2=.7*lambda2;end
      
   Xw=fft(fh);Xw2=abs(Xw).^2;
   Hw=(conj(Xw).*Yw)./(Xw2+lambda1*Dhw);
   hh=real(ifft(Hw));
   hh=hh(1:lh);
  
   hh=(1/sum(hh))*hh;
   iz=find(hh < 0);if(~isempty(iz)),hh(iz)=zeros(length(iz),1);end;
   hh=[hh;zeros(lg-lh,1)];

   Hw=fft(hh);Hw2=abs(Hw).^2;
   Xw=(conj(Hw).*Yw)./(Hw2+lambda2*Dfw2);
   fh=real(ifft(Xw)); 

   gh=H*fh;
  
   dg=g-gh;dhh=hh(2:lh)-hh(1:lh-1);
   J2=dg'*dg+lambda1*(dhh'*dhh)+lambda2*(fh'*fh);disp([k,J2])
   
   %J2=dg'*dg
   if(k>1),dJ=J2-J1;J1=J2;else,dJ=-1;J1=J2;end

   figure(1),clf;plot(f,'r'),hold on;plot(u1,'b'),plot(8*h1,'r'),plot(g,'b');axis([0 65 -.2 1.2]);xlabel('t')
   plot(fh,'+g');plot(fh,'g');th=text(32,.9,'f(t)');set(th,'Color','r'),th=text(37,.9,'fh(t)');set(th,'Color','g')
   plot(8*hh,'+c');plot(8*hh,'c');th1=text(5,.9,'h(t)');set(th1,'Color','r'),th2=text(10,.9,'hh(t)');set(th2,'Color','c')
   plot(gh,'+c');plot(gh,'c');th1=text(40,.7,'g(t)');set(th1,'Color','b'),th2=text(45,.7,'gh(t)');set(th2,'Color','c')
   hold off
   if(k==1),eval(['print -depsc Eps/xuhgh',num2str(nm1),'1.eps']);end 
  pause(1)
  end
  eval(['print -depsc Eps/xuhgh',num2str(nm1),'f.eps']);
 end
%---------------------------------------------------------------------------
 if (nm1 == 11)
  rep=menu('Satisfied ?','YES','NON');
  switch rep
  case 1
   disp('')
   disp('So, Thanks. Do not forget to mention it.')
   disp('If you want more, please contact Ali Mohammad-Djafari.')
   disp('Visit http://djafari.free.fr')
  case 2
   disp('')
   disp('So, Thank too.')
   disp('Can you do better ? Contact me. We can discuss!')
  end
  break
 end
%---------------------------------------------------------------------------
 if ((nm1>1)&(nm1<11))
  df=(f-fh);dg=g-gh;df=sum(df(:).*df(:))/sum(f(:).*f(:));
  dg=sum(dg(:).*dg(:))/sum(g(:).*g(:));
  disp('     df      dg ');disp([num2str(df), ' ',num2str(dg)]);
  %eval(['print -depsc Eps/xuhgh',num2str(nm1),'f.eps']);
 end

end

%---------------------------------------------------------------------------
function f=getnew(D,Text);
%GETNEW is equivalent to INPUT with default value
%f=getnew(D,Text);
%Test:
%f=1;f=getnew(f,'A value for x')
%Author: Ali Mohammad-Djafari
%
lD=length(D);
if(lD>1)
 TD='D: [';
 for k=1:lD
  TD=[TD,num2str(D(k)),';'];
 end
 TD(length(TD))=']';
 g=input([Text,' ? ',TD]);
else
 g=input([Text,' ? ','(D: ',num2str(D),') =>']);
end
[nx,ny]=size(g);
if (nx~=0)
 f=g; 
else
 f=D; 
end
disp(f)