function T = dttmtx(N,type,symv)
% DTTMTX Discrete trigonometric transform matrix (orthogonal)
%   T = dttmtx(N,type,symv)
%
%   This generates all 16 DTT transformation matrices.
%
%   N:    transform size (NxN)
%   type: dtt type. Here are all the possible types:
%         'dct1e','dct2e','dct3e','dct4e'
%         'dct1o','dct2o','dct3o','dct4o'
%         'dst1e','dst2e','dst3e','dst4e'
%         'dst1o','dst2o','dst3o','dst4o'
%   symv: if this flag is set then a symbolic matrix is produced
%         used for various proofs
%   T:    the transformation matrix
%
%   Notes: The implementation generates orthogonal DTT matrices.
%          as tabulated by Wang.  I am using ndgrid() here instead
%          outer products for readability.  Look at gdftmtx() for
%          the outer product implementation.
%
%   Example:
%   Transforming an framed input signal x (frame per column)
%   can be done as:
%   [N,M] = size(x);
%   T = dttmtx(N,'dct2e'); % Type II even
%   X = T*x;
%
%   Reference: Z. Wang and B. Hunt, The discrete W-transform
%              Appl. Math. Comput., 16 (1985), pp. 19-48.

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

% The default DTT type is DCT type II even
if nargin < 2; type = 'dct2e'; end
if nargin < 3; symv = 0;       end

% Make the indices n (time) and m (frequency) ...
[m,n]   = ndgrid(0:N-1);

% ... and the scaling factors
k      = ones(N,1);
k(1)   = 1/sqrt(2);
l      = ones(N,1);
l(end) = 1/sqrt(2);

% Make k,l diagonal matrices so that we can
% left / right multiply the kernel.  This is
% not computationaly efficient but it shows
% nicelly how the weights are applied.
k = diag(k);
l = diag(l);

% If the flag is setc convert stuff to
% symbolic and we are all set
if symv
    N = sym(N);
    m = sym(m); n = sym(n);
    k = sym(k); l = sym(l);
end

% Make the transformation matrix
switch lower(type)
    % ----------------------------------------------------
    % DCT I-IV even
    case {'dctie',  'dct1e'}
        N = N-1; k(end) = k(1);
        T = sqrt(2/N) *k *cos( m     .* n     *pi /N) *k;
    case {'dctiie', 'dct2e'}
        T = sqrt(2/N) *k *cos( m     .*(n+.5) *pi /N);
    case {'dctiiie','dct3e'}
        T = sqrt(2/N)    *cos((m+.5) .* n     *pi /N) *k;
    case {'dctive', 'dct4e'}
        T = sqrt(2/N)    *cos((m+.5) .*(n+.5) *pi /N);
    % ----------------------------------------------------
    % DCT I-IV odd
    case {'dctio',  'dct1o'}
        T = 2/sqrt(2*N-1) *k *cos( m     .* n     *2 *pi /(2*N-1)) *k;
    case {'dctiio', 'dct2o'}
        T = 2/sqrt(2*N-1) *k *cos( m     .*(n+.5) *2 *pi /(2*N-1)) *l;
    case {'dctiiio','dct3o'}
        T = 2/sqrt(2*N-1) *l *cos((m+.5) .* n     *2 *pi /(2*N-1)) *k;
    case {'dctivo', 'dct4o'}
        N = N+1;
        T = 2/sqrt(2*N-1)    *cos((m+.5) .*(n+.5) *2 *pi /(2*N-1));
    % ----------------------------------------------------
    % DST I-IV even
    case {'dstie',  'dst1e'}
        m = m+1; n = n+1; N = N+1;
        T = sqrt(2/N)    *sin( m     .* n     *pi /N);
    case {'dstiie', 'dst2e'}
        m = m+1; n = n+1; k(end) = k(1); k(1) = 1;
        T = sqrt(2/N) *k *sin( m     .*(n-.5) *pi /N);
    case {'dstiiie','dst3e'}
        m = m+1; n = n+1; k(end) = k(1); k(1) = 1;
        T = sqrt(2/N)    *sin((m-.5) .* n     *pi /N) *k;
    case {'dstive', 'dst4e'}
        T = sqrt(2/N)    *sin((m+.5) .*(n+.5) *pi /N);
    % ----------------------------------------------------
    % DST I-IV odd
    case {'dstio',  'dst1o'}
        m = m+1; n = n+1; N = N+1;
        T = 2/sqrt(2*N-1)    *sin( m     .* n     *2 *pi /(2*N-1));
    case {'dstiio', 'dst2o'}
        m = m+1; n = n+1; N = N+1;
        T = 2/sqrt(2*N-1)    *sin( m     .*(n-.5) *2 *pi /(2*N-1));
    case {'dstiiio','dst3o'}
        m = m+1; n = n+1; N = N+1;
        T = 2/sqrt(2*N-1)    *sin((m-.5) .* n     *2 *pi /(2*N-1));
    case {'dstivo', 'dst4o'}
        T = 2/sqrt(2*N-1) *l *sin((m+.5) .*(n+.5) *2 *pi /(2*N-1)) *l;
    % ----------------------------------------------------
    otherwise
        error('Unknown DTT type, try dct2e.  help dttmtx for a full list.')
end