function [f,R1,R2,err]=mpo(D,s,epsilon,k);

%         [f,R1,R2,err]=mpo(D,s,epsilon,k);
%
% Effectue le Matching-Pursuit Orthogonal (M.P.O.) d'un signal s de taille
% (1,N) suivant un dictionnaire D dont les atomes sont les D(j,:), j=1,...,
% m. Le dictionnaire D est une une matrice de taille (m,N). L'algorithme
% s'arrête au bout de k opérations ou dès que l'énergie E(r) du reste est
% telle que E(r) < epsilon*E(s), où E(s) désigne l'énergie du signal. La
% commande plot(R1) (respectivement plot(R2) donne (en abscisses) les
% numéros des atomes sélectionnés et (en ordonnées) la valeur de la partie
% réelle (respectivement imaginaire) du coefficient trouvé pour cet atome
% dans le resumé au cran k de l'algorithme. Les atomes sont supposés
% normalisés, d'énergie 1. On désigne par err l'erreur d'énergie relative.
% La variable de sortie f est une matrice de taille (m,N) dont les lignes
% successives, f(:,j), j=1,2,..., sont les approximations successives de s. 

[m,N]=size(D);
R1=zeros(1,m);
R2=zeros(1,m);
r=zeros(1,N);
E=s*s';
err=1;
err0=epsilon;
M=D*s';
M1=abs(M);
[Y,pmx]=max(M1);
a=M(pmx);
epmx=D(pmx,:)*D(pmx,:)';
newatom=D(pmx,:)/sqrt(epmx);
resume=(a/sqrt(epmx))*newatom;
coeff=[a/sqrt(epmx)];
natom=[pmx];
atoms=[newatom];
f=[resume];
invgram=[1];
D(pmx,:)=zeros(1,N);
for i=2:k
   if err0<err     
      reste=s-resume;
      M=D*reste';
      MM=abs(M);
      [Y,pmx]=max(MM);
      a=M(pmx);
      epmx=D(pmx,:)*D(pmx,:)';
      newatom=D(pmx,:)/sqrt(epmx);
      CC=newatom*atoms';
      C=conj(CC)';
      B=invgram*C;
      BB=conj(B');
      ed=1-(BB*CC');
      newcoeff=ed^(-1)*(reste*newatom');
      a1=coeff-newcoeff*BB;
      coeff=[a1,newcoeff];
      atoms=[atoms; newatom];
      resume=coeff*atoms;
      natom=[natom,pmx];
      D(pmx,:)=zeros(1,N);
      reste=s-resume;
      err=reste*reste'/E;
      f=[f;resume];
      beta=(1-conj(CC)*B)^(-1);
      invgram=[invgram + beta*B*B', -beta*B; -beta*B', beta]; 
   end  
end
k1=length(natom);
for j=1:k1, 
   l=fix(natom(j));
   R1(l)=real(coeff(j));
   R2(l)=imag (coeff(j));
end   

