function decmalp_b();
%            Deconvolution is an ill posed inverse problem
% f(t) input signal
% g(t) output signal 
% g(t)=h(t)*f(t)+e(t)
%
% Author: Ali Mohammad-Djafari
% Revision date: 10/1/2006
%
% 1. Load the data
% 2. No noise, everything works
% 3. With noise, no simple method work
% 4. Ad hoc methods
% 5. Wienner Filtering
% 6. Regularization Regularization (zero order)
% 7. Regularization Regularization (first order)
% 8. Quadratic Regularization: Iterative methods
% 9. Your method
% 0. Quit

clear all, close all, clc, help decmalp_b
figure(1);h=gcf;set(h,'Position',[500,300,300,300]);
figure(2);h=gcf;set(h,'Position',[700,300,300,300]);
figure(3);h=gcf;set(h,'Position',[900,300,300,300]);

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

 nm1=menu('Deconvolution is an ill posed inverse problem','Load the data', 'No noise: all methods work','With noise: no simple methods work','Ad hoc methods','Wienner filtering','Quadratic Regularization (zero order)','Qadratic Regularization (first order)','Qadratic Regularization: Iterative algorithm','Your method','Quit');
 switch nm1
 %---------------------------------------------------------------------------
 case 1
  disp('Load the data')
  load fhgb.mat

  figure(1),clf;plot(f,'r'),hold on;plot(f,'.r'),plot(lh*h1,'g'),plot(lh*h1,'.g'),plot(g,'b');plot(g,'.b');xlabel('t');axis([0 N 0 2]);
  th=text(N/2+10,1.2,'f(t)');set(th,'Color','r')
  th=text(5,.9,'h(t)');set(th,'Color','g')
  th=text(N/2+10,1.0,'g(t)=h(t)*f(t)');set(th,'Color','b')
  hold off
  print -depsc Eps/fhg.eps;

  figure(2),clf;plot(f,'r'),hold on;plot(f,'.r'),plot(lh*h1,'g'),plot(lh*h1,'.g'),plot(gb,'b');plot(gb,'.b');xlabel('t');axis([0 N 0 2]);
  th=text(N/2+10,1.2,'f(t)');set(th,'Color','r')
  th=text(5,.9,'h(t)');set(th,'Color','g')
  th=text(N/2+10,1.0,'gb(t)=h(t)*f(t)+e(t)');set(th,'Color','b')
  hold off
  print -depsc Eps/fhgb.eps;
 %---------------------------------------------------------------------------
 case 2
  load fhgb.mat
  disp('No noise : Any method may work')
  disp('Load data without noise')
  load fhgb.mat
  disp('Without noise even simple methods works. For example:')
  disp('  -- Use Inverse Filtering to do deconvolution.')

  Fh=G./(H+.001); 
  fh=real(ifft(Fh));
  Gh=H.*fft(fh);gh=real(ifft(Gh));
 %---------------------------------------------------------------------------
 case 3
  load fhgb.mat
  disp('With noise, no naive methods work')
  Fh=Gb./(H+.1); fh=real(ifft(Fh));
  Gh=H.*fft(fh);gh=real(ifft(Gh));
 %---------------------------------------------------------------------------
 case 4
  load fhgb.mat
  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));
  Fh=Gb./H;
  fh=real(ifft(Fh));
  Gh=H.*fft(fh);gh=real(ifft(Gh));
 %---------------------------------------------------------------------------
 case 5
  load fhgb.mat
  disp('Wienner Filtering') 
  r=20;r=getnew(r,'SNR in dB');  %r=10 20 30
  r1=10^(-r/20);

  Fh=Gb.*conj(H)./(abs(H).^2+r1); 
  fh=real(ifft(Fh));
  Gh=H.*fft(fh);gh=real(ifft(Gh));
 %---------------------------------------------------------------------------
 case 6
  load fhgb.mat
   disp('Quadratic Regularization (zero order): Try the values of lambda .1 1 10')
   lambda=.1;lambda=getnew(lambda,'Regularization coefficient lambda');
   H=toeplitz([h(1);zeros(lg-1,1)],h1);H=H';HtH=H'*H;
   DtD=eye(lg,lg);
   fh=(HtH+lambda*DtD)\(H'*gb);
   gh=conv(fh,h);gh=gh(1:lf);
 %---------------------------------------------------------------------------
 case 7
  load fhgb.mat
   disp('Quadratic Regularization (First order): Try the values of lambda .1 1 10')
   lambda=.1;lambda=getnew(lambda,'Regularization coefficient lambda');
   H=toeplitz([h(1);zeros(lg-1,1)],h1);H=H';HtH=H'*H;
   D=convmtx([1,-1],lg-1);DtD=D'*D;
   fh=(HtH+lambda*DtD)\(H'*gb);
   gh=conv(fh,h);gh=gh(1:lf);
 %---------------------------------------------------------------------------
 case 8
  load fhgb.mat
   disp('Quadratic Regularization: Iterative method')
   lambda=.1;lambda=getnew(lambda,'Regularization coefficient lambda');
   H=toeplitz([h(1);zeros(lg-1,1)],h1);H=H';Ht=H';HtH=H'*H;
   D=toeplitz([1;-1;zeros(lf-2,1)],[1;zeros(lf-1,1)]);Dt=D';DtD=Dt*D;

   itermax=50;fh=Ht*gb;alpha=1;
   for iter=1:itermax
    dg=gb-H*fh;df=D*fh;
    fh=fh+alpha*(Ht*dg-lambda*Dt*df);
    fh(fh<0)=0;
    figure(1),clf;plot(f,'r');hold on,plot(fh,'g');legend([' f';'fh']);axis([0 N -.2 2.2]);xlabel('t');hold off,pause(.1)
   end  
  gh=H*fh;
 %---------------------------------------------------------------------------
 case 9 % Your method
  load fhgb.mat
   disp('Your method')

  fh=my_name(g,H);

  gh=H*fh;
 end
 %---------------------------------------------------------------------------
 if((nm1>1)&(nm1<10))
  figure(1),clf;plot(f,'r');hold on,plot(fh,'g');legend([' f';'fh']);axis([0 N -.2 2.2]);xlabel('t');hold off
  eval(['print -depsc Eps/fh',num2str(nm1),'.eps']);
  figure(2),clf;plot(g,'b');hold on,plot(gh,'g');legend([' g';'gh']);axis([0 N -.2 2.2]);xlabel('t');hold off
  eval(['print -depsc Eps/gh',num2str(nm1),'.eps']);

  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)]);
 end
 if(nm1>9)
  break
 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)