% Best Approximation of a matrix by a frame multiplier
% function [TA uppsym] = ApprFramMult(T,D,Ds)
%
% usage:  TA=ApprFramMult(T,D)
%         [TA uppsym] = ApprFramMult(T,D)
%
% A matrix is best approximated (in the Hilbert-Schmidt sense) by a 
% frame multiplier. The elements of the frames are given in the synthesis
% matrices D and Ds columnswise. 
%  
% inputs : T .......... the matrix (m x n)
%          D .......... the elements of the analysis frame (columnwise) 
%                       (n x K)
%          Ds ......... the elements of the synthesis frame (columnwise) 
%                       (m x K)
%
% output: TA ......... the best apporximation of the matrix T with frame 
%                      multipliers using the frames in D and Ds
%         uppsym ..... symbol of best approximation
%
% test:   T = eye(2,2)
%         D = [0 1/sqrt(2) -1/sqrt(2); 1 -1/sqrt(2) -1/sqrt(2)]
%         [TA coeff] = ApprFramMult(T,D);
%         See testapprfarmmult.m and testapprGabmultKap1.m
%
% date: 10/03/2005 - 22/10/2007
%
% notes : D and Ds are just the synthesis operators of the respective frames. 
%        To be able to define frame multipliers they have to have the same  
%        number of elements. 
%
% Author:  XXL .. Peter Balazs, contact: a8927259@unet.univie.ac.at
%
% Literature : [1] P. Balazs; Irregular And Regular Gabor frame multipliers 
%                  with application to psychoacoustical masking 
%                  (Ph.D. thesis 2005)
%              [2] P. Balazs; Hilbert- Schmidt Operators and Frames -
%                  Classification, Best Approximation by Multipliers and 
%                  Algorithms; 
%                  International Journal of Wavelets, Multiresolution and
%                  Information Processing}, to appear, 
%                  http://arxiv.org/abs/math.FA/0611634
%
% See also: GMappir,low2uppIr
%
% Copyright :  (c) Acoustics Research Institute, Austrian Academy of
%              Science
%              http://www.kfs.oeaw.ac.at
function [TA,uppsym]=ApprFramMult(T,D,Ds)

if nargin < 2
    error('At leatst two inputs are needed; T and D');
end;
[N M] = size(T);
[Nd Kd] = size(D);

if N ~= Nd
    error('The number of rows in D and T have to be the same.');
    % as T:C^m -> C^n and D=( g_1 g_2 ...) with g_i \in C^n 
end;

if nargin < 3
    Ds = D;
else
    [Ns Ks] = size(Ds);
    if N ~= Ns
        error('The number of rows in Ds and T have to be the same.');
    end;
    if Kd ~= Ks
        error('The frames must have the same number of elements.');
    end;
end;

lowsym = zeros(Kd,1); %lower symbol
for i=1:Kd
%    d = D(:,i);
%    ds = Ds(:,i);
%    lowsym(i) = conj(d'*(T*ds));
    lowsym(i) = conj(Ds(:,i)'*(T*D(:,i)));
end;
% the more elegant 
% lowsym = diag(D'*T*D) 
% is slower, O(k(n^2+n^2)))
% see [Xxl]

% Gram-matrix in Hilbert-Schmidt sense
if nargin < 3
    Gram = abs((D'*D)).^2;
else
    Gram = (Ds'*Ds).*((D'*D).');
end;

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

% synthesis
TA = zeros(N,M);
for i = 1:Kd
%    d = D(:,i); 
%    ds = Ds(:,i);
%    P = d*ds';
    P = Ds(:,i)*D(:,i)';
    TA = TA + uppsym(i)*P;
end;
% found no faster or more elgant way






