function w=VSompJHH(x,S,Arr);
% VSompJHH  Orthogonal Matching Pursuit, Vector Selection algorithm
% as John H錵on Hus鴜 presented in lectures at HiS 1999.
% This algorithm orthogonalize each selected vector to the previous selected vectors.
%
% w=VSompJHH(x,S);
% w=VSompJHH(x,S,Arr);
%-----------------------------------------------------------------------------------
% arguments:
%   w    - the weights where at most S of the elements are non-zero, size Kx1
%   x    - the vector to approximate, size Nx1
%   S    - number of vectors to select or non-zero weights 
%          (do not mix with S in FindW function, which has a different meaning)
%   Arr  - acceptable energy of residual, (r=x-Fw; rr=r'*r;), if residual has lower
%          energy than Arr then function returns even if fewer than S vectors 
%          are selected, default is x'*x*1e-6
%-----------------------------------------------------------------------------------
% global variables:
%   F    - the normalized dictionary (or F matrix), size NxK
%          (this matrix is called U in some papers describing OMP)
%-----------------------------------------------------------------------------------

%-----------------------------------------------------------------------------------
% Copyright (c) 1999.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  05.11.1999  KS: function made
% Ver. 2.0  20.03.2000  KS: changed the arguments used 
% Ver. 2.1  18.12.2000  KS: use some global variables 
%-----------------------------------------------------------------------------------

global F
Mfile='VSompJHH';
[N,K]=size(F);
if nargin<2
   error([Mfile,': wrong number of arguments, see help.']);
end
if nargin<4
   Arr=x'*x*1e-6;
end

% start algorithm
epsilon=1e-8;       % a small number
r=x;                % the error so far
coeff=zeros(1,S);   % the coefficients
u=zeros(N,S);       % the orthogonalized vectors from F (up to S vectors)
n2u=zeros(1,S);     % norm squared of the u vectors
J=[];               % this will be the indexes of the selected vectors
k=0;                % number of vectors selected so far
while ((r'*r)>Arr)
   k=k+1;
   if k>S; break; end; 
   c=F'*r;   % the inner products
   [temp,i]=max(abs(c));i=i(1);
   J=[J,i];    % we select this vector from F, add index to J
   u(:,k)=F(:,i);
   for j=1:(k-1)    % orthogonalize u(:,k) vector to the vectors u(:,1:(k-1))
      u(:,k)=u(:,k)-((F(:,i)'*u(:,j))/n2u(j))*u(:,j);
   end
   n2u(k)=u(:,k)'*u(:,k);  % norm squared
   if n2u(k)<epsilon; n2u(k)=epsilon; end;  % we do not want to divide by zero
   coeff(k)=(u(:,k)'*r)/n2u(k);
   r=r-coeff(k)*u(:,k);
end

w=zeros(K,1);
if length(J)
   Fj=F(:,J);    % Fj is here only the selected vectors
   w(J)=((Fj'*Fj)\(Fj'))*x;
end

return;