[586fe9]: wavelets / wpbest.m Maximize Restore History

Download this file

wpbest.m    381 lines (329 with data), 10.8 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
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|.
%
% 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
if nargin<3
error('%s: Too few input parameters.',upper(mfilename));
end;
if ~isnumeric(J) || ~isscalar(J)
error('%s: "J" must be a scalar.',upper(mfilename));
end;
if(J<1 || rem(J,1)~=0)
error('%s: J must be a positive integer.',upper(mfilename));
end
if numel(nonzeros(size(f)>1))>1
error('%s: Function accepts only single channel inputs.',upper(mfilename));
end
definput.import = {'fwt'};
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.
wt = wfbtinit({w,J,'full'},'nat');
c = wpfbt(f,wt,'nat',flags.ext,'noscale');
% Nodes in the reverse BF order
treePath = nodesBForder(wt,'rev');
% Relationships between nodes
[pOutIdxs,chOutIdxs] = rangeWpBF(wt,'rev');
% Nodes to be removed
removeNodes = [];
% Energy normalization
% totalE = sum(cellfun(@(cEl) sum(cEl.^2),c));
% c = cellfun(@(cEl) cEl.^2/totalE,c,'UniformOutput',0);
if do_additive
% pre-calculate entropy of all subbands
cEnt = cellfun(@(cEl) wcostwrap(cEl,kv.cost,0),c);
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
cEnt = cell(numel(c),1);
for ii=1:length(pOutIdxs)-1
if isempty(cEnt{pOutIdxs(ii)})
pEnt = wcostwrap(c(pOutIdxs(ii)),kv.cost,0);
cEnt{pOutIdxs(ii)} = pEnt;
else
pEnt = cEnt(pOutIdxs(ii));
end
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
% Set parent entropy to the concatenation of the children entropy.
cEnt{pOutIdxs(ii)} = chEnt;
end
end
end
% Do tree prunning.
for ii=1:length(removeNodes)
wt = deleteSubtree(removeNodes(ii),wt);
end
end
% Finally do the analysis using the created best tree. with correct
% frequency bands order
[c,info] = wfbt(f,wt,flags.ext,'freq');
%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));