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