function f=mpo1(D,p,s,epsilon,k);

%           f=mpo1(D,p,s,epsilon,k);
%
% Effectue le matching-pursuit orthogonal (M.P.O.) d'un signal réel s de
% taille (1,N) suivant un dictionnaire D dont les atomes sont les D(j,:),
% j=1,..., m, et organisés en r sous-dictionnaires D1,...,Dr. La matrice D
% est 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. L'affichage des
% étiquettes des atomes retenus (ainsi que la valeur du coefficient trouvé
% au moment de la première détection de cet atome) se fait au fur et a
% mesure du deroulement de l'algorithme (les atomes sont tous supposés
% normalisés, d'énergie 1); les atomes sont classés au moment de
% l'affichage suivant le sous-dictionnaire auquels ils appartiennent. Faire
% [enter] après chaque affichage pour suivre cet affichage au fil de
% l'algorithme. La matrice 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. 

figure
[m,N]=size(D);
kk=fix(m/p);
R=zeros(1,kk); 
r=zeros(1,N);
E=s*s';
err=1;
err0=epsilon;
M=D*s';
M1=abs(M);
[Y,pmx]=max(M1);
pmx1=fix((pmx-1)./kk)+1;
pmx2=rem((pmx-1),kk)+1;
a=M(pmx);
epmx=D(pmx,:)*D(pmx,:)';
aa=a/sqrt(epmx);
newatom=D(pmx,:)/sqrt(epmx);
resume=(a/sqrt(epmx))*newatom;
coeff=[aa];
subplot(p,1,pmx1)
R(pmx2)=aa;
plot(R)
hold on
pause
R=zeros(1,kk);
atoms=[newatom];
f=[resume];
invgram=[1];
D(pmx,:)=zeros(1,N);
ll=1;
ll
for i=2:k
   if err0<err
      reste=s-resume;
      M=D*reste';
      MM=abs(M);
      [Y,pmx]=max(MM);
      pmx1=fix((pmx-1)./kk)+1;
      pmx2=rem(pmx-1,kk)+1;
      a=M(pmx);
      epmx=D(pmx,:)*D(pmx,:)';
      aa=a/sqrt(epmx);
      subplot(p,1,pmx1)
      R(pmx2)=aa;
      plot(R)
      hold on
      pause
      R=zeros(1,kk);
      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;
      D(pmx,:)=zeros(1,N);
      err=(s-resume)*(s-resume)'/E;
      f=[f;resume];
      beta=(1-conj(CC)*B)^(-1);
      invgram=[invgram + beta*B*B', -beta*B; -beta*B', beta];
      ll=ll+1;
      ll
   end  
end
   


