This is where navigation should be.

MULTIDGTREALMP - Matching Pursuit Decomposition with Multi-Gabor Dictionary

Usage

c = multidgtrealmp(f,dicts)
c = multidgtrealmp(f,dicts,errdb,maxit)
[c,frec,info] = multidgtrealmp(...)

Input parameters

f Input signal
dicts Dictionaries. Format {g1,a1,M1,g2,a2,M2,...}
errdb Target normalized approximation error in dB
maxit Maximum number of iterations.

Output parameters

c Sparse representation
frec Reconstructed signal
info Struct with additional output paramerets

Description

multidgtrealmp(f,{g1,a1,M1,g2,a2,M2,...,gW,aW,MW}) returns sparse representation of a signal in W Gabor dictionaries using the fast matching pursuit algorithm. gw is a Gabor window defined as in dgt and dgtreal, aw is a hop factor, Mw is the number of frequency channels. All aw and Mw must be divisible by min(a1,...,aW). The algorithm will try to reach -40 dB relative approximation error in at most numel(f) iterations. The function returns a cell-array with elements storing coefficients for individual Gabor systems such that they can be directly used in idgtreal.

multidgtrealmp(f,dicts,errdb,maxit) tries to reach normalized approximation error errdb dB in at most maxit iterations.

[c,frec,info] = multidgtrealmp(...) in addition returns the aproximation frec and a struct info with the following fields:

.iter Number of iterations done.
.atoms Number of atoms selected.
.relres Final approximation error. If resets are enabled, this is a vector additionally containing the approximation error at every reset.
.g Cell array of numeric windows used in the multi-dictionary
.a Array of hop factors for indivitual dictionaries
.M Array of numbers of channels for individual dictionaries
.synthetize Anonymous function which can be used to synthetize from the (modified) coefficients as frec = sum(info.synthetize(c),dim) where dim=2 if the input f was a column vector and dim=1 if it was a row vector.

The normalized approximation error is computed as errdb=20*log10(norm(f-frec)/norm(f)).

The function takes the following optional parameters at the end of the line of input arguments:

'kenrnthr',kt Kernel threshold. Must be in range ]0,1]. Default is 1e-4.
'timeinv' Use the time invariant phase convention. The default is 'freqinv'.
'pedanticsearch' Be pedantic about the energy of pairs of conjugated atoms in the selection step. Disabled by default.
'algorithm',alg Algorithm to use. Available: 'mp'(default),'selfprojmp','cyclicmp'
'reset' Reset the decomposition until stopping criteria are met. The reset means a complete synthesis and re-analysis. Available: 'noreset' (default), 'reset', 'quickreset' Option 'quickreset' applies a heuristic to stop the algorithm more quickly when the targer error is reached.

Reset conditions can be controlled by the following key-value pairs:

'resetit',it Reset the decomposition every it iteration. Must be a positive integer.
'reseterrdb',err Reset the decomposition when the error estimate in dB decreases below err. Must be negative. Default is 10*log10(1.1*kv.kernthr).

The computational routine is only available in C. Use ltfatmex to to compile it.

Algorithms

By default, the function uses the fast MP using approximate update in the coefficient domain as described in:

"Z. Prusa, N. Holighaus, P. Balazs: Fast Matching Pursuit with Multi-Gabor Dictionaries"

The kernel threshold limits the minimum approximation error which can be reached. For example the default threshold 1e-4 in general allows achieving at least -40 dB.

Note that this algorithm internally computes the error as errdb = 10*log10(err(k)/norm(f)^2), with err(0) = norm(f)^2 and err(k) = err(k-1) - abs(cselected(k))^2. Here, cselected(k) is the Gabor coefficient selected at the k-th MP step. This formula is exact only if the kernel threshold equals 0. Otherwise it serves as a cheap estimate. The exact error is only computed when this function is run with resets enabled and only in between resets.

Examples

The following example shows the decomposition in 3 dictionaries and plots contributions from the individual dictionaries and the residual.:

[f,fs] = gspi;
[c, frec, info] = multidgtrealmp(f,...
{'blackman',128,512,'blackman',512,2048,'blackman',2048,8192});
frecd = info.synthetize(c);
figure(1);
xvals = (0:numel(f)-1)/fs;
subplot(4,1,1); plot(xvals,frecd(:,1));ylim([-0.5,0.5]);
subplot(4,1,2); plot(xvals,frecd(:,2));ylim([-0.5,0.5]);
subplot(4,1,3); plot(xvals,frecd(:,3));ylim([-0.5,0.5]);
subplot(4,1,4); plot(xvals,f-frec);ylim([-0.5,0.5]); xlabel('Time (s)');

References:

S. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries. IEEE Trans. Signal Process., 41(12):3397--3415, 1993.

Z. Průša. Fast matching pursuit with multi-gabor dictionaries. Submitted., 2018. [ .pdf ]

L. Rebollo-Neira, M. Rozložník, and P. Sasmal. Analysis of a low memory implementation of the orthogonal matching pursuit greedy strategy. CoRR, abs/1609.00053v2, 2017.

B. L. Sturm and M. G. Christensen. Cyclic matching pursuits with multiscale time-frequency dictionaries. In Conf. Record of the 44th Asilomar Conference on Signals, Systems and Computers, pages 581--585, Nov 2010.