% Test the Approximation of Matrices by Gabor and Wavelet Mulitpliers
% Script-File
%
% This simple script tests the approximation of a random matrix 
% ('maximally unstructured system') by frame multipliers, using a Gabor 
% and wavelet frames. 
%
% date: 10/03/2005 - 22/10/2007 
%
% Author:  XXL .. Peter Balazs, contact: a8927259@unet.univie.ac.at
%
% Literature : [1] P. Balazs; Gabor frame multipliers with application to
%                  psychoacoustical masking (Ph.D. thesis) 
%              [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
%
% third party programs: 
%               gaussnk ... found in 'NuHAG-TB (BASIC!)' at
%                           http://www.univie.ac.at/nuhag-php/matlab/
%
% see also TestGabMulAppIrr2.m, GMAPPirr.m, testwavelet.m, waveletmat.m
%
% Copyright :  (c) Acoustics Research Institute, Austrian Academy of
%              Science, http://www.kfs.oeaw.ac.at

plotPlot = 1; % plot?

n = 32;
g = gaussnk(n); % analysis window
g = g*sqrt(0.7); % arbitrary nomalization factor to make frame bounds comparable
gamma = g; % synthesis window
a = 2; % Gabor parameter time shift
b = 4; % freq. shift
xpo = zeros(n,n);
xpo(1:b:n,1:a:n) = 1;

T = randc(n,n); % Target system

[TAI,COEFFI] = GMAPPir(T,xpo,g,g); % approximation by irregular Gabor algor.
G = gabbaspirr(g,xpo);
% S = G'*G;
disp('** Gabor system **');
if rank(G) ~= n
    error('Gabor atoms do not form frame!');
else
    [mGG nGG] = size(G);
    disp(sprintf('Redundancy is %g',mGG/n));
end
e = svd(G).^2;
% e = e(find(e>0));

disp(sprintf('Lower Frame Bound : A = %g Upper Frame Bound: B = %g',min(e),max(e)));

dilparam = 2; % wavelet parameter: dilation
shift = 1; % time shift
M = ceil(log(n)/log(dilparam))-1; % number of scales
a = [dilparam.^(-(M-2):-1) 1 dilparam.^(1:M+1)]

disp('** Wavelet system **');
WM = waveletmat(n,a,shift,1); % create wavelet synthesis matrix
if rank(WM) ~= n
    error('Wavelets do not form frame!');
end
[mWm nWm] = size(WM);
disp(sprintf('Redundancy is %g',mWm/n));
eWm = svd(WM);    
disp(sprintf(...  % continued on the next line for printout
' Lower Frame Bound A = %g Upper Frame Bound: B = %g',...
min(eWm),max(eWm)));
[appWM coeffWM] = ApprFramMult(T,WM.'); 
    % using standard Approximation by Frame Multiplier

TAC = bestapprcircl(T); % approximation by circular matrix = time-variant system.
if plotPlot == 1
 figure(1)
 subplot(2,2,1);
 surf(T);title('original');  axis([1 n 1 n -0.5 0.5]);
 subplot(2,2,2);surf(TAC); title('approximation by circulant matrix (time invariant system)');  axis([1 n 1 n -0.5 0.5]);
 subplot(2,2,3); surf(TAI); title('approximation by Gabor multipliers');  axis([1 n 1 n -0.5 0.5]);
 subplot(2,2,4); surf(appWM); title('approximation by wavelet multipliers');  axis([1 n 1 n -0.5 0.5]);
 
 figure(2);
 TS = fftshift(spread(T));
 TAIS = fftshift(spread(TAI));
 TACS = fftshift(spread(TAC));
 appWmsp = fftshift(spread(appWM));

 subplot(2,2,1);
 surf(TS);title('TF spread of original spread'); %axis([-ceil(n/2)-1 floor(n/2) -ceil(n/2)-1 floor(n/2) 0 8]);
   axis([1 n 1 n -4 4]);
   set(gca,'ytick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'yticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
   set(gca,'xtick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'xticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
 subplot(2,2,2);
 surf(TACS); title('TF spread of circular approximation');
   axis([1 n 1 n -4 4]);
   set(gca,'ytick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'yticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]); 
   set(gca,'xtick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'xticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
 subplot(2,2,3);
 surf(TAIS); title('TF spread of Gabor multiplier');
   axis([1 n 1 n -4 4]);
   set(gca,'ytick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'yticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
   set(gca,'xtick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'xticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
 subplot(2,2,4);
 surf(appWmsp); title('TF spread of wavelet multiplier');  
   axis([1 n 1 n -4 4]);
   set(gca,'ytick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'yticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
   set(gca,'xtick', [ceil(n/2-2*n/5) ceil(n/2-n/5) ceil(n/2) ceil(n/2+n/5) ceil(n/2+2*n/5)]);
   set(gca,'xticklabel', [ceil(-2*n/5) ceil(-n/5) 0 ceil(n/5) ceil(2*n/5)]);
end
 
DDG = diag(T)-diag(TAI); nDDG = norm(DDG);
DDW = diag(T)-diag(appWM); nDDW = norm(DDW);
disp(sprintf('Error of the approximation of diagonal: \n by Gabor multiplier (%f) and wavelet multiplier (%f). Factor: %f',nDDG,nDDW,nDDW/nDDG));



