matchPurs(x,W)

function [S,R,e,indx] = matchPurs(x,W)
% This function computes the projection of a given input vector or matrix
% onto a "dictionary" of other vectors or matrices using a matching pursuit
% algorithm.
%
% USAGES
% [S,R,e,indx] = matchPurs(x,W)  
% 
% INPUT
% x:	An Mx1 or MxN array. This array synthesized using dictionary 
%       elements from matrix W.
% W:    An MxN or MxNxP array of dictionary elements used to synthesize
%       input x. If x is an Mx1 vector, W must be an MxN matrix. If x is
%       an MxN array, W must be an MxNxP matrix containing the dictionary
%       elements.
%
% OUTPUT
% S:    The projection of each residual onto each dictionary element.
% R:    The residual x - sum(W,dim), where dim is 2 if x is a vector and
%       is 3 if x is a matrix.
% e:    The projection coefficients.
% indx  The index vector for the dictionary elements used in projections
%
% -------------------------------------------------------------------------
% TEST
% To ensure that the signal energy is preserved, the following relations
% should hold:
%
% Matrix Case:
% 
% (1) x = sum(S,3) + R;
% (2) norm(x,'fro') = sum(e.^2) + norm(R,'fro').^2
%
% Vector Case:
%
% (1) x = sum(S,2) + R;
% (2) norm(x).^2 = sum(e.^2) + norm(R).^2;


%% For Hilbert Spaces Whose Elements are Vectors
% This space has Euclidean inner product.  A vector is synthesized as:
% 
% x = ( x'*W(:,:,1) )*W(:,:,1) + (R1'*W(:,:,2) )*W(:,:,2) + ...
%     ( Rk'*W(:,:,k) )*W(:,:,k);
%
% Where Rk is the k_th residual.
% -------------------------------------------------------------------------

if( ~isvector(x) )
    
    [M,N,P] = size(W);
    S       = zeros(size(W));
    %intialize vectors
    R       = x; %the residual
    e       = zeros(P,1); %the energy vector
    ipvec   = zeros(P,1); %the dictionary indices
    indx    = ipvec; %the index for the best match
    
    for n = 1:P
        
        for k = 1:P
            
            ipvec(k)    = trace(R'*W(:,:,k))./norm( W(:,:,k), 'fro');
        
        end;
        
        [m, i]      = max(abs(ipvec));
        px          = ipvec(i)*W(:,:,i)./norm( W(:,:,i), 'fro');
        S(:,:,n)    = px; %keep each S
        R           = R - px;
        indx(n)     = i;
        W(:,:,i)    = nan(M,N);
        e(n)        = ipvec(i);
        
    end;
    
end;

% -------------------------------------------------------------------------
%% For Hilbert Spaces Whose Elements are Vectors
% This space has Euclidean inner product.  A vector is synthesized as:
% 
% x = ( x'*W(:,1) )*W(:,1) + (R1'*W(:,2) )*W(:,2) + ...
%     ( Rk'*W(:,k) )*W(:,k);
%
% Where Rk is the k_th residual.
% -------------------------------------------------------------------------

if( isvector(x) )
    
    [M,P]   = size(W);
    W       = normColumns(W);
    %intialize vectors
    R       = x; %the residual
    e       = zeros(P,1); %the energy vector
    ipvec   = zeros(P,1); %the dictionary indices
    indx    = ipvec; %the index for the best match
    
    for n = 1:P
        
        for k = 1:P
            
            ipvec(k)    = R'*W(:,k);
        
        end;
        
        [m, i]      = max(abs(ipvec));
        px          = ipvec(i)*W(:,i);
        S(:,n)      = px;
        R           = R - px;
        indx(n)     = i;
        W(:,i)      = nan(M,1);
        e(n)        = ipvec(i);
        
    end;
    
end;