function [c,info] = wpbest(f,w,J,varargin)
%WPBEST Best Tree selection
% Usage: c = wpbest(f,w,J,cost);
% [c,info] = wpbest(...);
%
% Input parameters:
% f : Input data.
% w : Wavelet Filterbank.
% J : Maximum depth of the tree.
%
% Output parameters:
% c : Coefficients stored in a cell-array.
% info : Transform parameters struct.
%
% `[c,info]=wpbest(f,w,J,cost)` selects the best sub-tree `info.wt` from
% the full tree with max. depth *J*, which minimizes the cost function.
%
% Only one-dimensional input *f* is accepted. The supported formats of
% the parameter *w* can be found in help for |fwt|. The format of the
% coefficients *c* and the `info` struct is the same as in |wfbt|.
%
% Please note that *w* should define orthonormal wavelet filters.
%
% First, the depth *J* wavelet packet decomposition is performed using |wpfbt|.
% Then the nodes are traversed in the breadth-first and bottom-up order
% and the value of the cost function of the node input and cost of the
% combined node outputs is compared. If the node input cost function value
% is less than the combined output cost, the current node and all
% possible descendant nodes are marked to be deleted, if not, the input is
% assigned the combined output cost. At the end, the marked nodes are
% removed and the resulting tree is considered to be a best basis (or
% near-best basis) in the chosen cost function sense.
%
% The `cost` parameter can be a cell array or an user-defined function handle.
% accepting a single column vector. The cell array should consist of a
% string, followed by a numerical arguments.
% The possible formats are listed in the following text.
%
% Additive costs:
% ---------------
%
% The additive cost *E* of a vector *x* is a real valued cost function
% such that:
%
% ..
% E(x) = sum E(x(k)),
% k
%
% .. math:: E(x)=\sum_{k} E(x(k))
%
% and $E(0)=0$. Given a collection of vectors $x_i$ being coefficients in
% orthonormal bases $B_i$, the best basis relative to *E* is the one for
% which the $E(x_i)$ is minimal.
%
% Additive cost functions allows using the fast best-basis search algorithm
% since the costs can be precomputed and combined cost of two vectors is
% just a sum of their costs.
%
% `{'shannon'}`
% A cost function derived from the Shannon entropy:
%
% ..
% E_sh(x) = -sum |x(k)|^2 log(|x(k)|^2),
% k:x(k)~=0
%
% .. math:: E_{sh}(x) = -\sum_{k:x(k)\neq 0} \left|x(k)\right|^2 \log \left|x(k)\right|^2
%
% `{'log'}`
% A logarithm of energy:
%
% ..
% E_log(x) = sum log(|x(k)|^2),
% k:x(k)~=0
%
% .. math:: E_{log}(x) = \sum_{k:x(k)\neq 0} \ln \left|x(k)\right|^2
%
% `{'lpnorm',p}`
% Concentration in $l^p$ norm:
%
% ..
% E_lp(x) = ( sum (|x(k)|^p) ),
% k
%
% .. math:: E_{lp}(x) = \left( \sum_{k} \left|x(k)\right|^p \right)
%
% `{'thre',th}`
% Number of coefficients above a threshold *th*.
%
%
% Non-additive costs:
% -------------------
%
% Cost function, which is not additive cost but which is used for the
% basis selection is called a non-additive cost. The resulting basis for
% which the cost is minimal is called near-best, because the non-additive
% cost cannot guarantee the selection of a best basis relative to the
% cost function.
%
% `{'wlpnorm',p}`
% The weak-$l^p$ norm cost function:
%
% ..
% E_wlp(x) = max k^{\frac{1}{p}}v_k(x),
%
% .. math:: E_{wlp}(x) = \max k^{1/p}v_k(x),
%
% where $0<p\leq 2$ and $v_k(x)$ denotes the *k*-th largest absolute value
% of *x*.
%
% `{'compn',p,f}`
% Compression number cost:
%
% ..
% E_cn(x) = arg min |w_k(x,p) - f|,
% k
%
% .. math:: E_{cn}(x) = \arg\min_{k} \left|w_k(x,p) - f\right|,
%
% where $0<p\leq 2$, $0<f<1$ and $w_k(u,p)$ denotes decreasingly sorted,
% powered, cumulateively summed and renormalized vector:
%
% .. k N
% w_k(x,p) = sum v_j^p(x) / sum v_j^p(x)
% j=1 j=1
%
% .. math:: w_k(x,p) = \frac{\sum_{j=1}^{k} v_j^p(x)}{\sum_{j=1}^{N} v_j^p(x)}
%
% where *v_k(x)* denotes the *k*-th largest absolute value of *x* and
% *N* is the number of elements of *x*.
%
% `{'compa',p}`
% Compression area cost:
%
% ..
% E_ca(x) = N - sum w_k(x,p),
% k
%
% .. math:: E_{ca}(x) = N - \sum_k w_k(x,p),
%
% where $0<p\leq 2$ and $w_k(u,p)$ and *N* as in the previous case.
%
% Examples:
% ---------
%
% A simple example of calling |wpbest| :::
%
% f = gspi;
% J = 8;
% [c,info] = wpbest(f,'sym10',J,'cost','shannon');
%
% % Use 2/3 of the space for the first plot, 1/3 for the second.
% subplot(3,3,[1 2 4 5 7 8]);
% plotwavelets(c,info,44100,90);
%
% subplot(3,3,[3 6 9]);
% N=cellfun(@numel,c); L=sum(N); a=L./N;
% plot(a,'o');
% xlim([1,numel(N)]);
% view(90,-90);
% xlabel('Channel no.');
% ylabel('Subsampling rate / samples');
%
% References: wickerhauser1991lectures tas94near-bestbasis
% AUTHOR: Zdenek Prusa
complainif_notenoughargs(nargin,3,'WPBEST');
complainif_notposint(J,'J','WPBEST');
if numel(nonzeros(size(f)>1))>1
error('%s: Function accepts only single channel inputs.',upper(mfilename));
end
definput.import = {'fwt','wfbtcommon'};
definput.flags.buildOrder = {'bottomup','topdown'};
definput.flags.bestWhat = {'tree','level'};
definput.keyvals.cost = {'shannon'};
[flags,kv]=ltfatarghelper({'cost'},definput,varargin);
if flags.do_topdown
error('%s: Flag %s not supported yet.',upper(mfilename),flags.buildOrder);
end
if flags.do_level
error('%s: Flag %s not supported yet.',upper(mfilename),flags.bestWhat);
end
if(~iscell(kv.cost))
kv.cost = {kv.cost};
end
% test if the chosen entropy measure is additive
do_additive = isAdditive(kv.cost);
if(flags.do_bottomup)
% Do full-tree Wavelet Packet decomposition beforehand and prune.
[c,info] = wpfbt(setnorm(f,'2'),{w,J,'full'},'nat',flags.ext,'intnoscale');
% Nodes in the reverse BF order
treePath = fliplr(nodeBForder(0,info.wt));
% Relationships between nodes
[pOutIdxs,chOutIdxs] = treeWpBFrange(info.wt);
% Nodes to be removed
removeNodes = [];
% Energy normalization
% totalE = sum(cellfun(@(cEl) sum(cEl.^2),c));
% c = cellfun(@(cEl) cEl.^2/totalE,c,'UniformOutput',0);
% pre-calculate entropy of all subbands
cEnt = cellfun(@(cEl) wcostwrap(cEl,kv.cost,0),c);
if do_additive
for ii=1:length(pOutIdxs)-1
pEnt = cEnt(pOutIdxs(ii));
chEnt = sum(cEnt(chOutIdxs{ii}));
if pEnt<=chEnt
removeNodes(end+1) = treePath(ii);
else
% Set parent entropy to the sum of the children entropy.
cEnt(pOutIdxs(ii)) = chEnt;
end
end
else
for ii=1:length(pOutIdxs)-1
pEnt = cEnt(pOutIdxs(ii));
chEnt = wcostwrap(c(chOutIdxs{ii}),kv.cost,0);
% Set parent entropy to value obtanied by concatenating child
% subbands.
if(pEnt<=chEnt)
removeNodes(end+1) = treePath(ii);
else
% Search the pOutIdxs(ii) in chOutIdxs
% There should be the only one..
foundId = cellfun(@(chEl)find(chEl==pOutIdxs(ii)),chOutIdxs,...
'UniformOutput',0);
chId = find(~cellfun(@isempty,foundId));
chOutIdxs{chId}(chOutIdxs{chId} == pOutIdxs(ii)) = [];
chOutIdxs{chId} = [chOutIdxs{chId},chOutIdxs{ii}];
end
end
end
% Do tree prunning.
for ii=1:length(removeNodes)
[info.wt] = nodeSubtreeDelete(removeNodes(ii),info.wt);
end
end
% Finally do the analysis using the created best tree. with correct
% frequency bands order
[c,info] = wfbt(f,info.wt,flags.ext,flags.forder);
%END WPBEST
function ad = isAdditive(entFunc)
x = 1:5;
%x = x./norm(x);
ent1 = wcostwrap(x',entFunc,0);
ent2 = sum(arrayfun(@(xEl) wcostwrap(xEl,entFunc,0),x));
ad = ent1==ent2;
function E = wcostwrap(x,fname,do_normalization)
%WENTWRAP Entropy functions wrapper
% Usage: E = wentwrap(x,fname,varargin)
%
% `E = wentwrap(x,fname,varargin)` passes given parameters further to the
% appropriate function.
%
if iscell(x)
x = cell2mat(x);
end
if do_normalization
x = x./norm(x);
end
if isa(fname,'function_handle')
E = fname(x);
return;
end
if numel(fname)>1
E = feval(sprintf('%s%s','wcost_',fname{1}),x,fname{2:end});
else
E = feval(sprintf('%s%s','wcost_',fname{1}),x);
end
%%%%%%%%%%%%%%%%%%
% ADDITIVE COSTS %
%%%%%%%%%%%%%%%%%%
function E = wcost_shannon(x)
% Cost function derived from the Shannon-Weaver entropy.
u = abs(x(x~=0)).^2;
E = -sum(u.*log(u));
function E = wcost_log(x)
% Logarithm of energy.
% It may be interpreted as the entropy of a Gauss-Markov process, composed
% of numel(x) uncorrelated Gaussian random variables of variences
% x(1),..,x(end). Minimizing the function finds the best approximation to
% the Karhuen-loéve basis for the process.
x = x(x~=0);
E = sum(log(x(:).^2));
function E = wcost_thre(x,th)
% Number of coefficients above a treshold.
% It gives a number of coefficients needed to transimt the signal to
% precision th.
E = numel(x(abs(x)>th));
function E = wcost_lpnorm(x,p)
% Concentration in l^p norm, p<2
% The smaller is the l^p norm of a signal with l^2 equal to 1, the more
% concentrated is its energy into a few coefficients.
assert(0<p&&p<2,'%s: Parameter p have to be in the range ]0,2[',upper(mfilename));
E = sum(abs(x(:)).^p);
%%%%%%%%%%%%%%%%%%%%%%
% NON-ADDITIVE COSTS %
%%%%%%%%%%%%%%%%%%%%%%
function E = wcost_wlpnorm(x,p)
% Weak Lp-norm
v = sort(abs(x(:)),'descend');
k = 1:length(v);
E = max((k.'.^(1/p)).*v);
function E = wcost_compa(x,p)
% Compresion Area Entropy, p<2
%
% 0<p<=2
assert(0<p&&p<=2,'Parameter p have to be in the range ]0,2]');
N = numel(x);
%u = x./norm(x);
v = sort(abs(x(:)),'descend');
wnom = cumsum(v.^p);
w = wnom./wnom(end);
E=N-sum(w);
function E = wcost_compn(x,p,f)
% Compresion Number entropy
% Represents a minimum number of coefficients containing 100*f% of the total
% p norm of x.
% 0<p<=2
assert(0<p&&p<=2,'Parameter p have to be in the range ]0,2]');
assert(0<f&&f<1,'Parameter f have to be in the range ]0,1[');
%u = x./norm(x);
v = sort(abs(x(:)),'descend');
wnom = cumsum(v.^p);
w = wnom./wnom(end);
[~,E]=min(abs(w-f));