function decmalpa();
%Impulse response identification and d�onvolution are ill-posed problems
%
% This demo shows the difficulties of:
% - Impulse response identification (x,y->h), 
% - Deconvolution (h,y->y) and 
% - Semi-Blind deconvolution (x1,y1->h and then h,y->x)
% - Blind deconvolution (y->h,x)
% First simple and naives methods are presented. 
% Then some regularization methods are applied.
%
% Author: Ali Mohammad-Djafari
% Revision date: 10/1/2007
% 
% Follow the menu.

% 1. Simulate experiences (generate h,x and y)
% 2. Apply simple and naives methods on no noisy data
% 3. Show that these naive methods cannot give satisfactory results 
%    when te data are noisy
% 4. Show that evoiding zero division methods are not satisfactory
% 5. Wienner Filtering methods 
% 6. Minimum Norm Least Square (MNLS) or zero order Regularization 
% 7. Regularization d'ordre un 
% 8. On estime h a partir de u et x par regularisation
% 9. Deconv aveugle: on estime h et x a partir de seul y
%10. Deconv aveugle: Mise en oeuvre par fft
% 0. Quit

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

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

 nm1=menu('Identification and Deconvolution are Ill-posed problemes','Simulate experiences', 'No noise case: Simple methods are OK','Noisy data: No satisfaction','Zero division evoiding methods','Wienner Filtering based methods','MNLS or Zero-order Regularization','Quadratic Regularization', 'Estimate h from u and x from y and he by regularization','Blind Deconvolution (y->x,h)','Blind Deconvolution: FFT based implementation','Quit');

 if (nm1 == 1)
  disp('1. Simulate experiences')
  nm2=menu('Simulate experiences','Example1','Example 2','Example 3','Example 4');
  switch nm2
  case 1
   t=[0:31]';
   x=[t;zeros(32,1)];x=(1/max(x))*x;lx=length(x);
   lh=32;a=[1;-1.6;.7];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);
   x=conv(x,h);x=x(3:lx+2);x=x/max(x);
  case 2
   x=[zeros(9,1);1;zeros(14,1);2;zeros(26,1);.5;zeros(12,1)];
   x=(1/max(x))*x;lx=length(x);
   lh=32;a=[1;-1.6;.7];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);
   x=conv(x,h);x=x(3:lx+2);x=x/max(x);lx=length(x);
  case 3
   x=[zeros(9,1);1;zeros(24,1);2;zeros(16,1);.5;zeros(12,1)];
   x=(1/max(x))*x;lx=length(x);
  case 4
   x=[zeros(5,1);ones(15,1);zeros(14,1);2*ones(14,1);zeros(16,1);];
   x=(1/max(x))*x;lx=length(x);
  end
  lh=32;a=[1;-1.2;.3];b=1;t=[1;zeros(lh-1,1)];h=filter(b,a,t);h=h/sum(h);
  y=conv(x,h);y=y(1:lx);ly=length(y);
  h1=[h;zeros(lx-lh,1)];lh1=length(h1);
  u=cumsum(h);lu=length(u);
  u1=cumsum(h1);lu1=length(u1);

  disp('   Save the data in xuhy.mat')
  save xuhy.mat x u h y h1 u1

  figure(1),clf;plot([x,u1,8*h1,y]);axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  h=text(33,.9,'x(t)');set(h,'Color','b')
  h=text(20,1.1,'u(t)');set(h,'Color','g')
  h=text(5,.9,'8*h(t)');set(h,'Color','r')
  h=text(40,.7,'y(t)');set(h,'Color','c')
  hold off
  %print -depsc eps/xuhy.eps;
  disp('   Look at Figure 1: h, u, x and y')
  disp('   Touch any key to continue'), pause(5)

  disp('   Adding noise to the data')
  yb=y+.01*randn(length(y),1);
  ub=u+.01*randn(length(u),1);
  u1b=u1+.01*randn(length(u1),1);
  disp('   Save the noisy data in xuhyb.mat')
  save xuhyb.mat x u h y h1 u1 yb ub u1b 

  figure(1),clf;plot([x,u1b,8*h1,yb]);axis([0 65 -.2 1.2]);xlabel('t')
  h=text(33,.9,'x(t)');set(h,'Color','b')
  h=text(20,1.1,'ub(t)');set(h,'Color','g')
  h=text(5,.9,'8*h(t)');set(h,'Color','r')
  h=text(40,.7,'yb(t)');set(h,'Color','c')
  hold on;
  %text(33,.9,'x(t)');text(20,1.1,'u(t)');text(5,.9,'8*h(t)');text(40,.7,'y(t)');
  hold off,pause(.1)
  %print -depsc eps/xuhyb.eps;
  disp('   Look at Figure 1: h, ub, x and yb')
 end

 if (nm1 == 2)
  disp('2. No noise case: Simple methods are OK')
  disp('   Load the data without noise')
  load xuhy.mat
  disp('   Simple methods works well: For example:')
  disp('    -- Derivating the step response u to obtain h')
  disp('    -- Use Inverse filtering to do deconvolution.')

  lx=length(x);ly=length(y);lh=length(h);lu=length(u);
  lh1=length(h1);lu1=length(u1);

  he=[0;u1(2:lu1)-u1(1:lu1-1)];
  H=fft(he); Y=fft(y); X=Y./(H+.1); 
  xe=real(ifft(X));
  ye=conv(xe,he);ye=ye(1:lx);

  figure(1),clf;plot([x,u1,8*h1,y]);axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  text(20,1.1,'u(t)');
  plot(xe,'+g');plot(xe,'g');text(32,.9,'x(t),  xe(t)');
  plot(8*he,'+b');plot(8*he,'b');text(5,.9,'h(t),  he(t)');
  plot(ye,'+c');plot(ye,'c');text(40,.7,'y(t), ye(t)');
  hold off,pause(.1)
  %eval(['print -depsc eps/xuhye',num2str(nm1),'.eps']);
 end

 if (nm1 > 2)
  disp('   Load the noisy data')
  load xuhyb.mat
  y=yb;
  u=ub;u1=u1b;
  lx=length(x);ly=length(y);lh=length(h);lu=length(u);
  lh1=length(h1);lu1=length(u1);
 end

 if (nm1 >= 3)&(nm1 <= 6)
  disp('   In this case, h is estimated from u by derivating.')
  he=[0;u1(2:lu1)-u1(1:lu1-1)];H=fft(he); 
 end

 if (nm1 == 3)
  disp('3. Noisy data : Simple methods do not work')
  Y=fft(y); X=Y./H; xe=real(ifft(X));
  ye=conv(xe,he);ye=ye(1:lx);
  
  figure(1),clf;plot([x,u1,8*h1,y]);axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  text(20,1.1,'u(t)');
  plot(xe,'+g');plot(xe,'g');text(32,.9,'x(t),  xe(t)');
  plot(8*he,'+b');plot(8*he,'b');text(5,.9,'h(t),  he(t)');
  plot(ye,'+c');plot(ye,'c');text(40,.7,'y(t),  ye(t)');
  hold off,pause(.1)
  %eval(['print -depsc eps/xuhye',num2str(nm1),'.eps']);
 end

 if (nm1 == 4)
  disp('4. Zero division evoiding methods are not satisfactoy') 
  eps1=.1;eps1=getnew(eps1,'   Threshold value');
  iz=find(abs(H)<eps1);H(iz)=eps1*ones(size(iz));
  Y=fft(y); 
  X=Y./H;
  xe=real(ifft(X));
  ye=conv(xe,he);ye=ye(1:lx);

  figure(1),clf;plot([x,y]);axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  plot(xe,'+g');plot(xe,'g');text(32,.9,'x(t),  xe(t)');
  plot(ye,'+c');plot(ye,'c');text(40,.7,'y(t), ye(t)');
  hold off,pause(.1)
  %eval(['print -depsc eps/xuhye',num2str(nm1),num2str(10*eps1),'.eps']);
 end

 if (nm1 == 5)
  disp('5. Wienner Filtering based methods') 
  r=20;r=getnew(r,'   SNR in dB');  %r=10 20 30
  r1=10^(-r/20);
  Y=fft(y); 
  X=Y.*conj(H)./(abs(H).^2+r1); 
  xe=real(ifft(X));
  ye=conv(xe,he);ye=ye(1:lx);

  figure(1),clf;plot([x,y]);axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  plot(xe,'+g');plot(xe,'g');text(32,.9,'x(t),  xe(t)');
  plot(ye,'+c');plot(ye,'c');text(40,.7,'y(t),  ye(t)');
  hold off,pause(.1)
  %eval(['print -depsc eps/xuhye',num2str(nm1),num2str(r/10),'.eps']);
 end

 if ((nm1 >= 6) & (nm1 < 8))
  he1=[0;u1(2:lu1)-u1(1:lu1-1)]; 
  H=toeplitz([zeros(ly,1)],he1);H=H'; 
  HtH=H'*H;
  if (nm1 == 6)
   disp('6. Zero-order Regularization: try different values of lambda .1 1 10')
   lambda=.1;lambda=getnew(lambda,'   Regul. parameter lambda');
   DtD=eye(ly,ly);
   xe=(HtH+lambda*DtD)\(H'*y);
  end
  if (nm1 == 7)
   disp('7. First-order Regularization: try different values of lambda .1 1 10')
   lambda=.1;lambda=getnew(lambda,'   Regul. parameter lambda');
   D=convmtx([1,-1],ly-1);DtD=D'*D;
   xe=(HtH+lambda*DtD)\(H'*y);
  end
  ye=conv(xe,he);ye=ye(1:lx);

  figure(1),clf;plot([x,y]);axis([0 65 -.2 1.2]);xlabel('t')
  text(40,.7,'y(t)');
  hold on;
  plot(xe,'+g');plot(xe,'g');text(32,.9,'x(t),  xe(t)');
  plot(ye,'+c');plot(ye,'c');text(40,.7,'y(t),  ye(t)');
  hold off,pause(.1)
  %eval(['print -depsc eps/xuhye',num2str(nm1),num2str(100*lambda),'.eps']);
 end
 
 if (nm1 == 8)
  disp('8. Semi Blind deconvolution')
  disp('   First estimate h using u and then use it for deconvolution')
  A=toeplitz([1;zeros(lu1-1,1)],ones(lu1,1));A=A'; 
  DtD=eye(ly,ly);
  lambda=1;lambda=getnew(lambda,'   Regul. parameter for h');
  he=(A'*A+lambda*DtD)\(A'*u1);
   
  disp('   Then estime x using he and y')
  H=toeplitz([he(1);zeros(ly-1,1)],he);H=H';HtH=H'*H; 
  lambda=.02;lambda=getnew(lambda,'   Regul. parameter for x');
  DtD=eye(ly,ly);
  xe=(HtH+lambda*DtD)\(H'*y);
  ye=conv(xe,he);ye=ye(1:lx);

  figure(1),clf;plot([x,u1,8*h1,y]);axis([0 65 -.2 1.2]);xlabel('t')
  hold on;
  text(20,1.1,'u(t)');text(40,.7,'y(t)');
  plot(xe,'+g');plot(xe,'g');text(32,.9,'x(t),  xe(t)');
  plot(8*he,'+b');plot(8*he,'b');text(5,.9,'h(t),  he(t)');
  plot(ye,'+c');plot(ye,'c');text(40,.7,'y(t),  ye(t)');
  hold off,pause(.1)
  %eval(['print -depsc eps/xuhye',num2str(nm1),'.eps']);
 end
 
 if (nm1 == 9)
  disp('9. Blind Deconvolution: y -> h, x')
  disp('   This is an iterative alg trying to optimize jointly')
  disp('   J(x,h)=|y-h*x|^2+lambda1*|D1 h|^2+lambda2*|D2 x|^2')
  disp('   0-Initialize xe and then iterate')
  disp('   1-Estimate he using xe and y')
  disp('   2-Estimate xe using he and y')
  disp('   At each iteration J is given:')
  lambda1=300;lambda1=getnew(lambda1,'   Regul. parameter for h');
  lambda2=1;lambda2=getnew(lambda2,'   Regul. parameter for x');
  xe=y;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([xe(1);zeros(lx-1,1)],xe);A=A'; 
   A=toeplitz([xe(1);zeros(lh-1,1)],xe);A=A'; 
   D=convmtx([1,-1],lh-1);DtD=D'*D;
   %D=eye(lh,lh);DtD=D;
   he=(A'*A+lambda1*DtD)\(A'*y); 
   he=(1/sum(he))*he;
   iz=find(he < 0);if(~isempty(iz)),he(iz)=zeros(length(iz),1);end;
   he1=[he;zeros(ly-lh,1)];
   dhe=D*he;
  
   %he(lh+1:lx)=zeros(lx-lh,1);
  
   H=toeplitz([he(1);zeros(ly-1,1)],he1);H=H'; 
   D=convmtx([1,-1],ly-1);DtD=D'*D;
   %DtD=eye(ly,ly);
   xe=(H'*H+lambda2*DtD)\(H'*y);
   dxe=D*xe;
   ye=H*xe;
  
   dy=y-ye;
   J2=dy'*dy+lambda1*(dhe'*dhe)+lambda2*(dxe'*dxe);
   disp(['    iter=',num2str(k),' J=',num2str(J2)]);

   %J2=dy'*dy
   if(k>1),dJ=J2-J1;J1=J2;else,dJ=-1;J1=J2;end

   plot([x,8*h1,y]);axis([0 65 -.2 1.2]);xlabel('t'),
   hold on
   plot(xe,'+g');plot(xe,'g');
   h=text(32,.9,'x(t)');set(h,'Color','b')
   h=text(36,.9,'xe(t)');set(h,'Color','g')
   plot(8*he,'+b');plot(8*he,'b');
   h=text(5,.9,'h(t)');set(h,'Color','g')
   h=text(9,.9,'he(t)');set(h,'Color','b')
   plot(ye,'+c');plot(ye,'c');
   h=text(40,.7,'y(t)');set(h,'Color','r')
   h=text(46,.7,'ye(t)');set(h,'Color','c')
   hold off,pause(1)
   %if(k==1),eval(['print -depsc eps/xuhye',num2str(nm1),'1.eps']);end
  end
  %eval(['print -depsc eps/xuhye',num2str(nm1),'f.eps']);
 end

 if (nm1 == 10)
  disp('10. Blind Deconvolution: y -> h, x')
  disp('   This is an iterative alg trying to optimize jointly')
  disp('   J(x,h)=|y-h*x|^2+lambda1*|D1 h|^2+lambda2*|D2 x|^2')
  disp('   0-Initialize xe and then iterate')
  disp('   1-Estimate he using xe and y')
  disp('   2-Estimate xe using he and y')
  disp('   At each iteration J is given:')
  disp('   Here steps 1 and 2 are implemented by FFT')
  lambda1=20;lambda1=getnew(lambda1,'   Regul. parameter for h');
  lambda2=.1;lambda2=getnew(lambda2,'   Regul. parameter for x');
  xe=y;lh=31;J1=0;J2=0;dJ=-1;
 
  Yw=fft(y);
  Dhw=fft([-1;1;zeros(ly-2,1)]);Dhw2=abs(Dhw).^2;
  Dfw=fft([-1;1;zeros(ly-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(xe);Xw2=abs(Xw).^2;
   Hw=(conj(Xw).*Yw)./(Xw2+lambda1*Dhw);
   he=real(ifft(Hw));
   he=he(1:lh);
  
   he=(1/sum(he))*he;
   iz=find(he < 0);if(~isempty(iz)),he(iz)=zeros(length(iz),1);end;
   he=[he;zeros(ly-lh,1)];

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

   ye=H*xe;
  
   dy=y-ye;dhe=he(2:lh)-he(1:lh-1);
   J2=dy'*dy+lambda1*(dhe'*dhe)+lambda2*(xe'*xe);
   disp(['    iter=',num2str(k),' J=',num2str(J2)]);
   
   %J2=dy'*dy
   if(k>1),dJ=J2-J1;J1=J2;else,dJ=-1;J1=J2;end

   plot([x,8*h1,y]);axis([0 65 -.2 1.2]);xlabel('t'),
   hold on
   plot(xe,'+g');plot(xe,'g');
   h=text(32,.9,'x(t)');set(h,'Color','b')
   h=text(36,.9,'xe(t)');set(h,'Color','g')
   plot(8*he,'+b');plot(8*he,'b');
   h=text(5,.9,'h(t)');set(h,'Color','g')
   h=text(9,.9,'he(t)');set(h,'Color','b')
   plot(ye,'+c');plot(ye,'c');
   h=text(40,.7,'y(t)');set(h,'Color','r')
   h=text(46,.7,'ye(t)');set(h,'Color','c')
   hold off,pause(1)
   %if(k==1),eval(['print -depsc eps/xuhye',num2str(nm1),'1.eps']);end 
  end
  %eval(['print -depsc eps/xuhye',num2str(nm1),'f.eps']);
 end

 if (nm1 == 11)
  disp('0. SATISFIED?')
  rep=menu('Satisfied ?','YES','NON');
  switch rep
  case 1
   disp('   So thanks. ')
   disp('   Do you want to know and to have more?, Contact Ali Mohammad-Djafari.')
  case 2
   disp('   Thanks too !')
   disp('   Can you do better? Contact me for discussions!')
  end
  disp('   Visit http://djafari.free.fr')
  break
 end

 if (1<nm1<11),
  %eval(['print -depsc eps/xuhye',num2str(nm1),'f.eps']);
 end

end

function x=getnew(D,Text);
%GETNEW is equivalent to INPUT with default value
%x=getnew(D,Text);
%Test:
%x=1;x=getnew(x,'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))=']';
 y=input([Text,' ? ',TD]);
else
 y=input([Text,' ? ','(D: ',num2str(D),') =>']);
end
[nx,ny]=size(y);
if (nx~=0)
 x=y; 
else
 x=D; 
end
disp(x)