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

%           [f,R1,R2,err]=mp(D,s,epsilon,k); 
%
% Effectue le Matching-Pursuit (M.P.) 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 matrice de taille (m,N). L'algorithme s'arrête
% après k itérations ou dès que l'énergie E(r) du reste est telle que l'on
% ait E(r) < epsilon * E(s), où E(s) désigne l'énergie du signal. 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.
% La commande plot(R1) (respectivement plot(R2)) donne les numéros des
% atomes sélectionnés (en abscisses) et les parties réelles (respectivement
% imaginaires) des coefficients dont ces atomes sont affectés (en
% ordonnées), lorsque les atomes sont normalisés et d'énergie égale à 1. La
% variable err désigne l'erreur relative commise après k itérations en
% termes d'énergie. 

[m,N]=size(D);
f=[];
resume=zeros(1,N);
E=s*s';
err=1;
err0=epsilon;
R1=zeros(1,m);
R2=zeros(1,m);
for i=1:k
   if err0<err
      M=D*s';
      M1=abs(M);
      [Y,pmx]=max(M1); 
      a=M(pmx);
      ed=D(pmx,:)*D(pmx,:)';
      R1(pmx)=R1(pmx)+ (real(M(pmx))/sqrt(ed));
      R2(pmx)=R2(pmx)+ (imag(M(pmx))/sqrt(ed));
      resume=resume+(a/ed)*D(pmx,:);
      s=s-(a/ed)*D(pmx,:);
      err=s*s'/E;
      f=[f;resume];    
   end  
end