% Calculation of cross-Gram Matrix for Gabor systems 
%  HSGramMatrXXL.M   XXL,  13.05.2004  
%
% function   hsgm = HsGramMatrXXL(xpo,g,gamma,show,full)
%
%  determines either the full Hilbert Schmidt Gram Matrix of the
%  rank one operators  f -> <f, \pi(\lambda) \gamma> \pi(\lambda')
%  g (full = 1) or only the entries important for a Gabor 
%  multiplier (full = 0)
%
%  inputs : xpo   .... (0/1) square matrix of points in the TF plane
%                      row index - frequency, column index -time 
%                      (following the convention for normal spectograms)
%           g     .... synthesis window row vector (1 x m), 
%                      Default: Gauss
%           gamma .... analysis window row vector (1 x m'), 
%                      Default: g
%           show  .... show graphics (slower) (*/0), Default: 1 
%                      (if input ~= 0, show = 1!)
%           full  .... calculate the HS Gram Matrix of all  
%                      possible TF projections 
%                      or only 
%                      the entries which are used in Gabor 
%                      multipliers, where the tensor product only
%                      depends on one TF point (1/*), 
%                      Default: 0 (if input ~= 1, full = 0!)
%
%  output: the k^2 x k^2 HS gram matrix, where k is the number of 
%          non zero entries in xpo. (full = 1)
%          the k x k part of HS gram matrix, important for Gabor 
%          multiplier (otherwise)
%
%  usage:  HsGramMatrXXL(xpo,g,gamma,show,full)
%
%  notes: The entries of this HS Gram Matrix are (full == 1)
%         < g_k \otimes gamma*_j , g_i \otimes gamma*_l >_HS =
%       = < g_k \otimes gamma*_j , g_i \otimes gamma*_l >_L2 =
%       = < g_k , g_i > < gamma_j, gamma_l >* =
%       = STFT_gamma(gamma)_(k,i) * STFT*_g(g)_(j,l)
%      
%       (full == 0) => k = j, i = l
%
% See also: GMAPPirr, teststft
%
% Ref.: [xxl]... P. Balazs, "Regular And Irregular Gabor Multipliers With
%                Application To Psychoacoustic Masking", PhD thesis, Univ.
%                of Vienna, 2005
% 
% complexity: O(n^4) (full == 1) or O(n^2) 
%             (so can be time consuming)
%
%  XXL .. Peter Balazs, contact: a8927259@unet.univie.ac.at
%
% notes: 
% The matrices A,B are (up to a phase factor) the "Cross"-Gram  
% Matrices of the TF atoms g_j, g_l resp. gamma_k, gamma_i.
%
% For full == 1, the complexity is about O(n^4), most of the time
% is used by the function "kron" if n is big enough. This function
% is a built in function, so should be fairly optimized.
% Otherwise most of the time is used up by the assignement of the
% matrcies.
% For full == 0, the complexity is O(n^2).
% For n=144, k = 0,5 % of N (full == 0) this functions needs about
% 12 sec. (PIII 937)
% For n=144, k = 1,5 % of N (full == 0) this functions needs about
% 101 sec.
%
% Copyright :  (c) NUHAG, Faculty of Math., University of Vienna, AUSTRIA
%              http://nuhag.mat.univie.ac.at/
%              (c) Acoustics Research Institute, Austrian Academy of
%              Science
%              http://www.kfs.oeaw.ac.at


function   hsgm = HsGramMatrXXL(xpo,g,gamma,show,full)
% to test how much time is used by which part, use profiling:
% profile on -detail builtin;
if nargin < 1
  error('Function HSGramMatrXXL has to get at least one parameter : xpo');
end
[N M] = size(xpo);
if M ~= N
    error('XPO has to be a square matrix');
end;
[xpx xpy] = find(xpo > 0);
% xpx ... frequencies
% xpy ... time
k = size(xpx,1);
if k == 0 
  error('The matrix xpo is containing only zeros');
end
if nargin < 2  
  g = gaussnk(N);
  disp('using default gaussian analysis window');
else
  % Maybe create a subfunction check(g);
  [n m] = size(g);
  if m < n 
    g = g.';
    disp('transposing g');
    [n m] = size(g);
  end
  if n ~= 1 
    error('Input g should be a vector')
  end  
  if m > N
    error('The size of g is bigger than the (spectral) size of xpo')
  elseif m < N    
    gaga = zeros(1,N - size(g)) % zeropadding
    gamma = [gamma;gaga]
  end
end
    
if nargin < 3
  gamma = g; 
  disp('using default synthesis window: gamma = g');
else
  [n m] = size(gamma);
  if m < n 
    gamma = gamma.';
     disp('transposing gamma');
     [n m] = size(gamma);
  end
  if n ~= 1 
    error('Input gamma should be a vector')
  end  
  if m > N
    error('The size of g is bigger than the (spectral) size of xpo')
  elseif m < N    
    gg = zeros(1,N - m) % zeropadding
    g = [g;gg]
  end
end

if nargin < 4
  show = 0;
elseif show ~= 0
  show = 0;
end

if nargin < 5
  full = 0;
elseif full ~= 1
  full = 0;
end
    
gst = stft(g,g,1,1); % full stft
gast = stft(gamma,gamma,1,1); % full stft

if show == 1
  figure(1);
  subplot(2,1,1); 
  imgc(gst); 
  title('Wigner-like distribution of analysis atom');
  subplot(2,1,2); 
  imgc(gast); 
  title('Wigner-like distribution of synthesis atom'); 
end
A = zeros(k);
B = zeros(k);

% try to programm this part using MATLAB matrix functions
for ii = 1:k
  for jj = 1:k
     A(ii,jj) = gst(mod(xpx(jj)-xpx(ii),N)+1,mod(xpy(jj)-xpy(ii),N)+1);  
     % row index of STFT are time values, see teststft.m
     % using periodic extension
     B(ii,jj) = conj(gast(mod(xpx(jj)-xpx(ii),N)+1,mod(xpy(jj)-xpy(ii),...
         N)+1)); 
     % conjugate for a 'real' tensor product !!!!!
    end
end

if full == 1 
  hsgm = kron(A,B);
else
  hsgm = A.*B;
end

if show == 1
  figure(2);
  SURF(abs(hsgm));
end
% Profiling end part:
% profile report;
% figure(3);
% profile plot;
% profile off;

