1function [tgrad,fgrad,s,c]=filterbankphasegrad(f,g,a,L,minlvl)
2%-*- texinfo -*-
3%@deftypefn {Function} filterbankphasegrad
4%@verbatim
5%FILTERBANKPHASEGRAD   Phase gradient of a filterbank representation
6%   Usage:  [tgrad,fgrad,s,c] = filterbankphasegrad(f,g,a,L,minlvl);
7%           [tgrad,fgrad,s,c] = filterbankphasegrad(f,g,a,L);
8%           [tgrad,fgrad,s,c] = filterbankphasegrad(f,g,a,minlvl);
9%           [tgrad,fgrad,s,c] = filterbankphasegrad(f,g,a);
10%           [tgrad,fgrad,s] = filterbankphasegrad(...)
11%           [tgrad,fgrad]  = filterbankphasegrad(...)
12%
13%   Input parameters:
14%      f     : Signal to be analyzed.
15%      g     : Cell array of filters
16%      a     : Vector of time steps.
17%      L     : Signal length (optional).
18%      minlvl: Regularization parameter (optional, required < 1).
19%   Output parameters:
20%      tgrad : Instantaneous frequency relative to original position.
21%      fgrad : The negative of the local group delay.
22%      cs    : Filterbank spectrogram.
23%      c     : Filterbank coefficients.
24%
25%   [tgrad,fgrad,s,c] = FILTERBANKPHASEGRAD(f,g,a,L) computes the
26%   relative instantaneous frequency tgrad and the negative of the group
27%   delay fgrad of the filterbank spectrogram s obtained from the
28%   signal f and filterbank parameters g and a.
29%   Both tgrad and fgrad are specified relative to the original
30%   coefficient position entirely similar to GABPHASEGRAD.
31%   fgrad is given in samples, while tgrad is given in normalised
32%   frequencies such that the absolute frequencies are in the range of ]-1,1].
33%
34%   This routine uses the equivalence of the filterbank coefficients in
35%   each channel with coefficients obtained from an STFT obtained with a
36%   certain window (possibly different for every channel). As a consequence
37%   of this equivalence, the formulas derived in the reference apply.
38%
39%
40%   References:
41%     F. Auger and P. Flandrin. Improving the readability of time-frequency
42%     and time-scale representations by the reassignment method. IEEE Trans.
43%     Signal Process., 43(5):1068--1089, 1995.
44%
45%     N. Holighaus, Z. Průša, and P. L. Soendergaard. Reassignment and
46%     synchrosqueezing for general time-frequency filter banks, subsampling
47%     and processing. Signal Processing, 125:1--8, 2016. [1]http ]
48%
49%     References
50%
51%     1. http://www.sciencedirect.com/science/article/pii/S0165168416000141
52%
53%@end verbatim
54%@strong{Url}: @url{http://ltfat.github.io/doc/filterbank/filterbankphasegrad.html}
55%@seealso{gabphasegrad}
56%@end deftypefn
57
58% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
59% This file is part of LTFAT version 2.3.1
60%
61% This program is free software: you can redistribute it and/or modify
62% it under the terms of the GNU General Public License as published by
63% the Free Software Foundation, either version 3 of the License, or
64% (at your option) any later version.
65%
66% This program is distributed in the hope that it will be useful,
67% but WITHOUT ANY WARRANTY; without even the implied warranty of
68% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
69% GNU General Public License for more details.
70%
71% You should have received a copy of the GNU General Public License
72% along with this program.  If not, see <http://www.gnu.org/licenses/>.
73
74%   AUTHOR : Nicki Holighaus.
75
76complainif_notenoughargs(nargin,3,'FILTERBANKPHASEGRAD');
77
78% Reshape input signal
79[f,~,W]=comp_sigreshape_pre(f,'FILTERBANKPHASEGRAD',0);
80Ls = size(f,1);
81
82if W>1
83    error('%s: Only one-channel signals supported.',upper(mfilename));
84end
85
86if nargin < 5
87    if nargin < 4
88        L = filterbanklength(Ls,a);
89        minlvl = eps;
90    else
91        if ~(isscalar(L) && isnumeric(L) ) && L>0
92            error('%s: Fourth argument shoud be a positive number.',...
93                  upper(mfilename));
94        end
95        if L >= 1
96            minlvl = eps;
97        else
98            minlvl = L;
99        end;
100    end
101end;
102
103complainif_notposint(L,'L','FILTERBANKPHASEGRAD');
104
105Luser = filterbanklength(L,a);
106if Luser~=L
107    error(['%s: Incorrect transform length L=%i specified. ', ...
108           'Next valid length is L=%i. See the help of ',...
109           'FILTERBANKLENGTH for the requirements.'],...
110           upper(mfilename),L,Luser);
111end
112
113
114% Unify format of coefficients
115[g,asan]=filterbankwin(g,a,L,'normal');
116
117% Precompute filters
118[gh, gd, g] = comp_phasegradfilters(g, asan, L);
119
120f=postpad(f,L);
121
122c=comp_filterbank(f,g,asan);
123% Compute filterbank coefficients with frequency weighted window
124ch=comp_filterbank(f,gh,asan);
125% Compute filterbank coefficients with time weighted window
126cd=comp_filterbank(f,gd,asan);
127
128% Run the computation
129[tgrad,fgrad,s] = comp_filterbankphasegrad(c,ch,cd,L,minlvl);
130
131