function F = filt_decomp(H,k,tol)

% filt_decomp - singular value decomposition of filter
% ----------------------------------------------------
%
% F = filt_decomp(H,k,tol)
%
% Input:
% ------
%  H - filter mask
%  k - rank of decomposition (def: 0, economy size full decomposition)
%  tol - rank computation tolerance (def: same as in 'rank')
%
% Output:
% -------
%  F - filter decomposition structure

%--------------------------------
% Author: Harold Figueroa
%--------------------------------
% $Date: 2004-12-02 15:58:47 -0500 (Thu, 02 Dec 2004) $
% $Revision: 132 $
%--------------------------------

%-------------------------------
% HANDLE INPUT
%-------------------------------

%--
% set order of decomposition
%--

if (nargin < 2)
	k = 0;
end

%--
% check order of decomposition
%--

[m,n] = size(H);

if (k > 0)
	if (k > min(m,n))
		k = 0;
	end
end

%-------------------------------
% COMPUTE DECOMPOSITION
%-------------------------------

%--
% compute full or partial singular value decomposition
%--

if (k == 0)
	[U,S,V] = svd(H,'econ');
else
	[U,S,V] = svds(H,k);
end

S = diag(S)';

%-------------------------------
% OUTPUT FILTER
%-------------------------------

F = filter_create;

%--
% set filter mask and decomposition fields
%--

F.H = H;

F.S = S; F.X = V; F.Y = U;

%--
% set tolerance based rank
%--

if (nargin < 3)
	F = filt_rank(F);
else
	F = filt_rank(F,tol);
end

%--
% approximation and computation fields
%--

% NOTE: the approximation perspective uses norm of remaining singular vectors

if (k == 0)
	E = 0;
else
	E = norm(H - (U * diag(S) * V'),'fro');
end

E = sqrt(cumsum([E, fliplr(S.^2)])); 

E(end) = []; E = fliplr(E);

F.error = E;

% NOTE: the theoretical speed is off due to overhead

d = length(S);

F.speed = (2*m*n - 1) ./ (((1:d) .* ((2*m - 1) + (2*n - 1))) + (0:(d - 1)));
