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
|
function sr=gabreassignadjust(s,pderivs,a,varargin)
%-*- texinfo -*-
%@deftypefn {Function} gabreassignadjust
%@verbatim
%GABREASSIGNADJUST Adjustable reassignment of a time-frequency distribution
% Usage: sr = gabreassignadjust(s,pderivs,a,mu);
%
% GABREASSIGNADJUST(s,pderivs,a,mu) reassigns the values of the positive
% time-frequency distribution s using first and second order phase
% derivatives given by pderivs and parameter mu*>0.
% The lattice is determined by the time shift a and the number of
% channels deduced from the size of s.
%
% pderivs is a cell array of phase derivatives which can be obtained
% as follows:
%
% pderivs = gabphasederiv({'t','f','tt','ff','tf'},...,'relative');
%
% Please see help of GABPHASEDERIV for description of the missing
% parameters.
%
% gabreassign(s,pderivs,a,mu,despeckle) works as above, but some
% coeficients are removed prior to the reassignment process. More
% precisely a mixed phase derivative pderivs{5} is used to determine
% which coefficients m,n belong to sinusoidal components (such that
% abs(1+pderivs{5}(m,n)) is close to zero) and to impulsive
% components (such that abs(pderivs{5}(m,n)) is close to zero).
% Parameter despeckle determines a threshold on the previous quantities
% such that coefficients with higher associated values are set to zeros.
%
% Algorithm
% ---------
%
% The routine uses the adjustable reassignment presented in the
% references.
%
% Examples:
% ---------
%
% The following example demonstrates how to manually create a
% reassigned spectrogram.:
%
% % Compute the phase derivatives
% a=4; M=100;
% [pderivs, c] = gabphasederiv({'t','f','tt','ff','tf'},'dgt',bat,'gauss',a,M,'relative');
%
% % Reassignemt parameter
% mu = 0.1;
% % Perform the actual reassignment
% sr = gabreassignadjust(abs(c).^2,pderivs,a,mu);
%
% % Display it using plotdgt
% plotdgt(sr,a,143000,50);
%
%
% References:
% F. Auger, E. Chassande-Mottin, and P. Flandrin. On phase-magnitude
% relationships in the short-time fourier transform. Signal Processing
% Letters, IEEE, 19(5):267--270, May 2012.
%
% F. Auger, E. Chassande-Mottin, and P. Flandrin. Making reassignment
% adjustable: The Levenberg-Marquardt approach. In Acoustics, Speech and
% Signal Processing (ICASSP), 2012 IEEE International Conference on,
% pages 3889--3892, March 2012.
%
% Z. Průša. STFT and DGT phase conventions and phase derivatives
% interpretation. Technical report, Acoustics Research Institute,
% Austrian Academy of Sciences, 2015.
%
%@end verbatim
%@strong{Url}: @url{http://ltfat.github.io/doc/gabor/gabreassignadjust.html}
%@seealso{gabphasederiv, gabreassign}
%@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: Peter L. Soendergaard, 2008; Zdeněk Průša 2015
thisname = upper(mfilename);
complainif_notenoughargs(nargin,3,thisname);
complainif_notposint(a,'a',thisname);
definput.keyvals.mu=0;
definput.keyvals.despeckle=0;
[~,~,mu,despeckle] = ltfatarghelper({'mu','despeckle'},definput,varargin);
if ~(isscalar(mu) && mu>=0)
error('%s: mu must be a real positive number.',thisname);
end
if ~(isscalar(despeckle) && despeckle>=0)
error('%s: despeckle must be a real positive number.',thisname);
end
[M,N,W] = size(s);
if W>1
error(['%s: c must be 2D matrix.'],thisname);
end
if ~(iscell(pderivs) && numel(pderivs) == 5)
error(['%s: pderiv must be a cell array of phase derivatives in ',...
'the following order t,f,tt,ff,tf.'],thisname);
end
% Basic checks
if any(cellfun(@(el) isempty(el) || ~isnumeric(el),{s,pderivs{:}}))
error(['%s: s and elements of the cell array pderivs must be ',...
'non-empty and numeric.'],upper(mfilename));
end
% Check if argument sizes are consistent
sizes = cellfun(@size,pderivs,'UniformOutput',0);
if ~isequal(size(s),sizes{:})
error(['%s: s and all elements of the cell array pderivs must ',...
'have the same size.'], upper(mfilename));
end
% Check if any argument is not real
if any(cellfun(@(el) ~isreal(el),{s,pderivs{:}}))
error('%s: s and all elements of the cell array pderivs must be real.',...
upper(mfilename));
end
if any(s<0)
error('%s: s must contain positive numbers only.',...
upper(mfilename));
end
[tgrad,fgrad,ttgrad,ffgrad,tfgrad] = deal(pderivs{:});
if despeckle~=0
% Removes coefficients which are neither sinusoidal component or
% impulse component based on the mixed derivative.
% How reassigned time position changes over time
thatdt = -tfgrad;
% How reassigned frequency position changes along frequency
ohatdo = 1+tfgrad;
% Only coefficients with any of the previous lower than despeckle is
% kept.
s(~(abs(ohatdo)<despeckle | abs(thatdt)<despeckle)) = 0;
end
% Construct the inverses explicitly
%
% |trelpos| = |A1 A2|^-1|B1|
% |frelpos| = |A3 A4| |B2|
%
% det(A)*|trelpos| = | A4 -A2|*|B1|
% |frelpos| = |-A3 A1 |B2|
B1 = fgrad(:);
B2 = tgrad(:);
A1 = tfgrad(:) + 1 + mu;
A2 = -ffgrad(:);
A3 = -ttgrad(:);
A4 = -tfgrad(:) + mu;
dets = (A1.*A4-A2.*A3);
oneoverdets=1./dets;
% Remove nearly singular matrices
% The coefficients will not be reassigned
oneoverdets(abs(dets)<1e-10) = 0;
trelpos = oneoverdets.*( A4.*B1 - A2.*B2);
frelpos = oneoverdets.*(-A3.*B1 + A1.*B2);
% frelpos is derived from tgrad and
% trelpos is derived from fgrad
sr=comp_gabreassign(s,reshape(frelpos,M,N),reshape(trelpos,M,N),a);
|
