function test_dtt
% Here we test all the DTTs by checking the inversion rules
% and the orthogonality

% ------- test_dtt.m ---------------------------------------
% Marios Athineos, marios@ee.columbia.edu
% http://www.ee.columbia.edu/~marios/
% Copyright (c) 2005 by Columbia University.
% All rights reserved.
% ----------------------------------------------------------

% We will test one odd and one even length signal
N = 10:19;

% Tolerance (-14 fails for size 72x72 and up)
% Checked with 10^-13 up to 1024x1024 and it's fine
tol = 10^-13;

% I make a cell array of all the types to check them in a loop
t = { 'dct1e','dct2e','dct3e','dct4e', ...
      'dct1o','dct2o','dct3o','dct4o', ...
	  'dst1e','dst2e','dst3e','dst4e', ...
	  'dst1o','dst2o','dst3o','dst4o' }';

  tic
% Lets test
for n = N
    fprintf('For N = %d, time %e\n',n,toc);
    tic
    % In order to get the kernel we use eye().  This
    % way we also test both dtt() as well as dttmtx()
    E = eye(n);
    
    for I = 1:length(t)
        % This is for inversion rules for orhogonal
        R = E - dtt(E,t{I},1) * idtt(E,t{I},1);
        if ~isempty(find(abs(R) > tol))
            fprintf('Inverse of %s (    orth) is wrong, I = %d\n',t{I},I);
        end
        
        % This is for orthogonality
        R = E - dtt(E,t{I},1) * dtt(E,t{I},1)';
        if ~isempty(find(abs(R) > tol))
            fprintf('%s is not orthogonal, I = %d\n',t{I},I);
        end

        % This is for inversion rules for non orthogonal
        R = E - dtt(E,t{I},0) * idtt(E,t{I},0);
        if ~isempty(find(abs(R) > tol))
            fprintf('Inverse of %s (non orth) is wrong, I = %d\n',t{I},I);
        end
    end
end