1function [xo]=groupthresh(xi,lambda,varargin)
2%-*- texinfo -*-
3%@deftypefn {Function} groupthresh
4%@verbatim
5%GROUPTHRESH   Group thresholding
6%   Usage:  xo=groupthresh(xi,lambda);
7%
8%   GROUPTHRESH(x,lambda) performs group thresholding on x, with
9%   threshold lambda.  x must be a two-dimensional array, the first
10%   dimension labelling groups, and the second one labelling members. This
11%   means that the groups are the row vectors of the input (the vectors
12%   along the 2nd dimension).
13%
14%   Several types of grouping behaviour are available:
15%
16%    GROUPTHRESH(x,lambda,'group') shrinks all coefficients within a given
17%     group according to the value of the l^2 norm of the group in
18%     comparison to the threshold lambda. This is the default.
19%
20%    GROUPTHRESH(x,lambda,'elite') shrinks all coefficients within a
21%     given group according to the value of the l^1 norm of the
22%     group in comparison to the threshold value lambda.
23%
24%   GROUPTHRESH(x,lambda,dim) chooses groups along dimension
25%   dim. The default value is dim=2.
26%
27%   GROUPTHRESH accepts all the flags of THRESH to choose the
28%   thresholding type within each group and the output type (full / sparse
29%   matrix). Please see the help of THRESH for the available
30%   options. Default is to use soft thresholding and full matrix output.
31%
32%
33%
34%   References:
35%     M. Kowalski. Sparse regression using mixed norms. Appl. Comput. Harmon.
36%     Anal., 27(3):303--324, 2009.
37%
38%     M. Kowalski and B. Torresani. Sparsity and persistence: mixed norms
39%     provide simple signal models with dependent coefficients. Signal, Image
40%     and Video Processing, 3(3):251--264, 2009.
41%
42%     G. Yu, S. Mallat, and E. Bacry. Audio Denoising by Time-Frequency Block
43%     Thresholding. IEEE Trans. Signal Process., 56(5):1830--1839, 2008.
44%
45%@end verbatim
46%@strong{Url}: @url{http://ltfat.github.io/doc/sigproc/groupthresh.html}
47%@seealso{thresh, demo_audioshrink}
48%@end deftypefn
49
50% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
51% This file is part of LTFAT version 2.3.1
52%
53% This program is free software: you can redistribute it and/or modify
54% it under the terms of the GNU General Public License as published by
55% the Free Software Foundation, either version 3 of the License, or
56% (at your option) any later version.
57%
58% This program is distributed in the hope that it will be useful,
59% but WITHOUT ANY WARRANTY; without even the implied warranty of
60% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
61% GNU General Public License for more details.
62%
63% You should have received a copy of the GNU General Public License
64% along with this program.  If not, see <http://www.gnu.org/licenses/>.
65
66%   AUTHOR : Kai Siedenburg, Bruno Torresani.
67%   REFERENCE: OK
68
69
70if nargin<2
71  error('Too few input parameters.');k
72end;
73
74if (prod(size(lambda))~=1 || ~isnumeric(lambda))
75  error('lambda must be a scalar.');
76end;
77
78% Define initial value for flags and key/value pairs.
79definput.import={'thresh','groupthresh'};
80definput.importdefaults={'soft'};
81definput.keyvals.dim=2;
82
83[flags,keyvals,dim]=ltfatarghelper({'dim'},definput,varargin);
84
85% kv.dim (the time or frequency selector) is handled by assert_sigreshape_pre
86[xi,L,NbMembers,NbGroups,dim,permutedsize,order]=assert_sigreshape_pre(xi,[],dim,'GROUPTHRESH');
87
88if flags.do_sparse
89  xo = sparse(size(xi));
90else
91  xo = zeros(size(xi));
92end;
93
94if flags.do_group
95
96  groupnorm = sqrt(sum(abs(xi).^2));
97  w = thresh(groupnorm, lambda, flags.iofun,flags.outclass)./groupnorm;
98
99  % Clean w for NaN. NaN appears if the input has a group with norm
100  % exactly 0.
101  w(isnan(w)) = 0;
102
103  xo = bsxfun(@times,xi,w);
104
105end
106
107if flags.do_elite
108  for ii=1:NbGroups,
109    y = sort(abs(xi(:,ii)),'descend');
110    rhs = cumsum(y);
111    rhs = rhs .* lambda ./ (1 + lambda * (1:NbMembers)');
112    M_ii = find(diff(sign(y-rhs)));
113    if (M_ii~=0)
114      tau_ii = lambda * norm(y(1:M_ii),1)/(1+lambda*M_ii);
115    else
116      tau_ii = 0;
117    end
118
119    % FIXME: The following line does not work for sparse matrices.
120    xo(:,ii) = thresh(xi(:,ii),tau_ii,flags.iofun,flags.outclass);
121  end
122end;
123
124xo=assert_sigreshape_post(xo,dim,permutedsize,order);
125
126
127