1function [sym,lowb,upb]=gabmulappr(T,p2,p3,p4,p5);
2%-*- texinfo -*-
3%@deftypefn {Function} gabmulappr
4%@verbatim
5%GABMULAPPR  Best Approximation by a Gabor multiplier
6%   Usage:  sym=gabmulappr(T,a,M);
7%           sym=gabmulappr(T,g,a,M);
8%           sym=gabmulappr(T,ga,gs,a,M);
9%           [sym,lowb,upb]=gabmulappr( ... );
10%
11%   Input parameters:
12%         T     : matrix to be approximated
13%         g     : analysis/synthesis window
14%         ga    : analysis window
15%         gs    : synthesis window
16%         a     : Length of time shift.
17%         M     : Number of channels.
18%
19%   Output parameters:
20%         sym   : symbol
21%
22%   sym=GABMULAPPR(T,g,a,M) calculates the best approximation of the given
23%   matrix T in the Frobenius norm by a Gabor multiplier determined by the
24%   symbol sym over the rectangular time-frequency lattice determined by
25%   a and M.  The window g will be used for both analysis and
26%   synthesis.
27%
28%   GABMULAPPR(T,a,M) does the same using an optimally concentrated, tight
29%   Gaussian as window function.
30%
31%   GABMULAPPR(T,gs,ga,a) does the same using the window ga for analysis
32%   and gs for synthesis.
33%
34%   [sym,lowb,upb]=GABMULAPPR(...) additionally returns the lower and
35%   upper Riesz bounds of the rank one operators, the projections resulting
36%   from the tensor products of the analysis and synthesis frames.
37%
38%
39%
40%   References:
41%     M. Doerfler and B. Torresani. Representation of operators in the
42%     time-frequency domain and generalized Gabor multipliers. J. Fourier
43%     Anal. Appl., 16(2):261--293, April 2010.
44%
45%     P. Balazs. Hilbert-Schmidt operators and frames - classification, best
46%     approximation by multipliers and algorithms. International Journal of
47%     Wavelets, Multiresolution and Information Processing, 6:315 -- 330,
48%     2008.
49%
50%     P. Balazs. Basic definition and properties of Bessel multipliers.
51%     Journal of Mathematical Analysis and Applications, 325(1):571--585,
52%     January 2007.
53%
54%     H. G. Feichtinger, M. Hampejs, and G. Kracher. Approximation of
55%     matrices by Gabor multipliers. IEEE Signal Procesing Letters,
56%     11(11):883--886, 2004.
57%
58%@end verbatim
59%@strong{Url}: @url{http://ltfat.github.io/doc/operators/gabmulappr.html}
60%@seealso{framemulappr, demo_gabmulappr}
61%@end deftypefn
62
63% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
64% This file is part of LTFAT version 2.3.1
65%
66% This program is free software: you can redistribute it and/or modify
67% it under the terms of the GNU General Public License as published by
68% the Free Software Foundation, either version 3 of the License, or
69% (at your option) any later version.
70%
71% This program is distributed in the hope that it will be useful,
72% but WITHOUT ANY WARRANTY; without even the implied warranty of
73% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
74% GNU General Public License for more details.
75%
76% You should have received a copy of the GNU General Public License
77% along with this program.  If not, see <http://www.gnu.org/licenses/>.
78
79% AUTHOR    : Monika Doeerfler
80% REFERENCE : REF_GABMULAPPR
81% TESTING   : TEST_GABMULAPPR
82
83complainif_argnonotinrange(nargin,3,5,mfilename);
84
85L=size(T,1);
86
87if size(T,2)~=L
88  error('T must be square.');
89end;
90
91if nargin==3
92  % Usage: sym=gabmulappr(T,a,M);
93  a=p2;
94  M=p3;
95  ga=gabtight(a,M,L);
96  gs=ga;
97end;
98
99if nargin==4
100  % Usage: sym=gabmulappr(T,g,a,M);
101  ga=p2;
102  gs=p2;
103  a=p3;
104  M=p4;
105end;
106
107if nargin==5
108  % Usage: sym=gabmulappr(T,ga,gm,a,M);
109  ga=p2;
110  gs=p3;
111  a=p4;
112  M=p5;
113end;
114
115if size(ga,2)>1
116  if size(ga,1)>1
117    error('Input g/ga must be a vector');
118  else
119    % ga was a row vector.
120    ga=ga(:);
121  end;
122end;
123
124if size(gs,2)>1
125  if size(gs,1)>1
126    error('Input g/gs must be a vector');
127  else
128    % gs was a row vector.
129    gs=gs(:);
130  end;
131end;
132
133N=L/a;
134b=L/M;
135
136Vg=dgt(gs,ga,1,L);
137
138s=spreadfun(T);
139
140A=zeros(N,M);
141V=zeros(N,M);
142for k=0:b-1
143  for l=0:a-1
144    A = A+ s(l*N+1:(l+1)*N,k*M+1:k*M+M).*conj(Vg(l*N+1:(l+1)*N,k*M+1:k*M+M));
145    V = V+abs(Vg(l*N+1:(l+1)*N,k*M+1:k*M+M)).^2;
146  end;
147end;
148
149if nargout>1
150  lowb = min(V(:));
151  upb  = max(V(:));
152end;
153
154SF1=A./V;
155
156SF=zeros(N,M);
157jjmod=mod(-M:-1,M)+1;
158iimod=mod(-N:-1,N)+1;
159SF=SF1(iimod,jjmod);
160
161sym=b*dsft(SF)*sqrt(M)/sqrt(N);
162
163