% Approximation of Matrices by irregular Gabor multiplier 
% function [TAI,COEFF]=gmappirr(T,xpo,g,gamma,show)
%
% The best approximation of a matrix by a Gabor multiplier with 
% analysis window g and synthesis window gamma on the (irregular) grid xpo 
% is calculated.
%  
% inputs : T ......... the matrix
%          g ......... Gabor analysis window, row vector
%          gamma ..... Gabor synthesis window, row vector
%          xpo ....... the (possibly) irregular grid, row index = frequency
%          show ...... flag if error should be calculated 
%
% output: TAI ....... the best approximation by Gabor multipliers
%         coeff ..... the lower symbol
%
% usage:  [TAI,COEFF] = Blo2WalXXL(A,n)
%
% test:   use TestGabMulAppIrr.m
%
% last change: 21/09/2007
%
% third party programs: 
%               gaussnk ... found in 'NuHAG-TB (BASIC!)' at
%                           http://www.univie.ac.at/nuhag-php/matlab/
%
% notes: The convention of the NuHAG tollboxes using row vectors and 
%        matrix multiplication from the right is used here !!
%
% see also: ApprFramMult.m, HSGramMatrxxl.m, testhsgrammatr.m
%
% Author:  XXL .. Peter Balazs, contact: a8927259@unet.univie.ac.at
%
% Literature : [Xxl] P. Balazs; Gabor frame multipliers with application to
%                    psychoacoustical masking 
%                    (Ph.D. thesis, in preparation) 
%              [FHK] H. G. Feichtinger, M. Hampjes, G. Kracher; 
%                    Approximation of Matrices by Gabor Multipliers, 
%                    IEEE Signal Procesing Letters Vol. 11, No. 11 (2004)
%              [Doe] M. Doerfler; Gabor Analysis for a Class of Signals  
%                    called Music, PhD thesis Univ. Wien (2002)
%
%
% 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 [TAI,uppsym]=gmappirr(T,xpo,g,gamma,show)

if nargin < 2
  error('Function GMAPPIr has to get at least two parameters : T,xpo');
end
[N M] = size(xpo);
xxx = find(xpo.' > 0);
k = size(xxx,1);
if k == 0 
  error('The matrix xpo is containing only zeros');
end

[TN TM] = size(T);
if TN ~= TM 
    error('At the moment only square matrices can be used.');
elseif TN ~= N
    error('Matrix size does not fit grid size.');
end;

if nargin < 3  
  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 | m < M    
    error('The length of g is too small for xpo.');
  end
end
    
if nargin < 4
  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 | m < M    
    error('The length of g is too small for xpo.');
  end
end

GBIa = gabbaspIrr(g,xpo);
GBIs = gabbaspIrr(gamma,xpo);

lowsym = zeros(k,1); %lower symbol
for ii=1:k
    lowsym(ii) = (GBIa(ii,:)*T)*(GBIs(ii,:)');
end;
% Gram-matrix in Hilbert-Schmidt sense
Gram = HSGramMatrXXL(xpo,g,gamma);

% upper symbol:
uppsym = pinv(Gram)*lowsym;

% synthesis
% TAI = zeros(N,M);
% for ii = 1:k
%     P = GBIa(ii,:)'*GBIs(ii,:);
%     TAI = TAI + uppsym(ii)*P;
% end;

TAI = (GBIa'*diag(uppsym(:)))*GBIs;

if nargin < 5 | show ~= 1
   show=0;
end
if show==1 
    disp('Fehler:');
    norm(T-TA,'fro')
end;

 
