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
|
function [xo,N]=thresh(xi,lambda,varargin);
%-*- texinfo -*-
%@deftypefn {Function} thresh
%@verbatim
%THRESH Coefficient thresholding
% Usage: x=thresh(x,lambda,...);
% [x,N]=thresh(x,lambda,...);
%
% THRESH(x,lambda) will perform hard thresholding on x, i.e. all
% elements with absolute value less than scalar lambda will be set to zero.
%
% THRESH(x,lambda,'soft') will perform soft thresholding on x,
% i.e. lambda will be subtracted from the absolute value of every element
% of x.
%
% The lambda parameter can also be a vector with number of elements
% equal to numel(xi) or it can be a numeric array of the same shape
% as xi. lambda is then applied element-wise and in a column major
% order if lambda is a vector.
%
% [x,N]=THRESH(x,lambda) additionally returns a number N specifying
% how many numbers where kept.
%
% THRESH takes the following flags at the end of the line of input
% arguments:
%
% 'hard' Perform hard thresholding. This is the default.
%
% 'wiener' Perform empirical Wiener shrinkage. This is in between
% soft and hard thresholding.
%
% 'soft' Perform soft thresholding.
%
% 'full' Returns the output as a full matrix. This is the default.
%
% 'sparse' Returns the output as a sparse matrix.
%
% The function wTHRESH in the Matlab Wavelet toolbox implements some of
% the same functionality.
%
% The following code produces a plot to demonstrate the difference
% between hard and soft thresholding for a simple linear input:
%
% t=linspace(-4,4,100);
% plot(t,thresh(t,1,'soft'),'r',...
% t,thresh(t,1,'hard'),'.b',...
% t,thresh(t,1,'wiener'),'--g');
% legend('Soft thresh.','Hard thresh.','Wiener thresh.','Location','NorthWest');
%
%
% References:
% S. Ghael, A. Sayeed, and R. Baraniuk. Improved wavelet denoising via
% empirical Wiener filtering. In Proceedings of SPIE, volume 3169, pages
% 389--399. San Diego, CA, 1997.
%
% J. Lim and A. Oppenheim. Enhancement and bandwidth compression of noisy
% speech. Proceedings of the IEEE, 67(12):1586--1604, 1979.
%
%@end verbatim
%@strong{Url}: @url{http://ltfat.github.io/doc/sigproc/thresh.html}
%@seealso{largestr, largestn}
%@end deftypefn
% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
% This file is part of LTFAT version 2.3.1
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR : Kai Siedenburg, Bruno Torresani and Peter L. Soendergaard.
% TESTING: OK
% REFERENCE: OK
complainif_notenoughargs(nargin,2,'THRESH');
is_sameshape = ndims(lambda)==ndims(xi) && all(size(lambda)==size(xi));
if ~isnumeric(lambda) || ...
~isscalar(lambda) && ... % lambda is not scalar
numel(lambda)~=numel(xi) && ... % lambda does not have the same number of elements
~(is_sameshape) % lambda does not have the same shape
error(['%s: lambda must be a scalar, a vector with ',...
'numel(lambda)==numel(xi) or whatever shape xi has such that ',...
'all(size(lambda)==size(xi))'],upper(mfilename));
end;
% Define initial value for flags and key/value pairs.
definput.import={'thresh'};
[flags,keyvals]=ltfatarghelper({},definput,varargin);
if flags.do_sparse
if ndims(xi)>2
error(['%s: Sparse output is only supported for 1D/2D input. This ',...
'is a limitation of Matlab/Octave.'],upper(mfilename));
end;
if ~isa(xi,'double')
error(['%s: Input is not double prec. data array and sparse output can,'...
'be double precision data type only. This is a ',...
'Matlab/Octave limitation.'],upper(mfilename));
end
end;
% Reshape lambda if it is a vector
if ~is_sameshape && ~isscalar(lambda)
lambda = reshape(lambda,size(xi));
end
if flags.do_sparse
xo=sparse(size(xi,1),size(xi,2));
if flags.do_hard
if isscalar(lambda)
% Create a significance map pointing to the non-zero elements.
signifmap=find(abs(xi)>=lambda);
else
signifmap=abs(xi)>=lambda;
end
xo(signifmap)=xi(signifmap);
else
if isscalar(lambda)
% Create a significance map pointing to the non-zero elements.
signifmap=find(abs(xi)>lambda);
else
signifmap=abs(xi)>lambda;
end
if flags.do_wiener
if isscalar(lambda)
xo(signifmap) = 1 - (lambda./abs(xi(signifmap))).^2;
else
xo(signifmap) = 1 - (lambda(signifmap)./abs(xi(signifmap))).^2;
end
xo(signifmap) = xi(signifmap).*xo(signifmap);
end;
if flags.do_soft
if isscalar(lambda)
% xo(signifmap)=xi(signifmap) - sign(xi(signifmap))*lambda;
xo(signifmap)=(abs(xi(signifmap)) - lambda) .* ...
exp(i*angle(xi(signifmap)));
else
xo(signifmap)=(abs(xi(signifmap)) - lambda(signifmap)) .* ...
exp(i*angle(xi(signifmap)));
end
% The line above produces very small imaginary values when the input
% is real-valued. The next line fixes this
if isreal(xi)
xo=real(xo);
end;
end;
end
if nargout==2
N=numel(signifmap);
end;
else
% Dense case
xo=zeros(size(xi),assert_classname(xi));
% Create a mask with a value of 1 for non-zero elements. For full
% matrices, this is faster than the significance map.
if flags.do_hard
if nargout==2
mask=abs(xi)>=lambda;
N=sum(mask(:));
xo=xi.*mask;
else
xo=xi.*(abs(xi)>=lambda);
end;
end;
if flags.do_soft
% In the following lines, the +0 is significant: It turns
% -0 into +0, oh! the joy of numerics.
if nargout==2
xa=abs(xi)-lambda;
mask=xa>=0;
xo=(mask.*xa+0).*sign(xi);
N=sum(mask(:))-sum(xa(:)==0);
else
xa=abs(xi)-lambda;
xo=((xa>=0).*xa+0).*sign(xi);
end;
end;
if flags.do_wiener
xa = lambda./abs(xi);
xa(isinf(xa)) = 0;
xa = 1 - xa.^2;
if nargout==2
mask = xa>0;
xo = xi.*xa.*mask;
N = sum(mask(:));
else
xo = xi.*xa.*(xa>0);
end
end;
end;
|
