function T = dttmtxno(N,type,frwd,symv)
% DTTMTX Discrete trigonometric transform matrix (non-orthogonal)
%   T = dttmtxno(N,type,frwd,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'
%   frwd: if set the forward transform is calculated otherwise
%         the inverse.  this is necessary cause on non-orthogonality
%   symv: if this flag is set then a symbolic matrix is produced
%         used for various proofs
%   T:    the transformation matrix
%
%   Notes: The implementation generates non orthogonal DTT matrices.
%          as defined by Martucci
%
%   Example:
%   Transforming an framed input signal x (frame per column)
%   can be done as:
%   [N,M] = size(x);
%   T = dttmtxno(N,'dct2e'); % Type II even
%   X = T*x;
%
%   Reference: S.A. Martucci, "Symmetric convolution and the discrete
%              sine and cosine transforms", IEEE Trans. Sig. Processing
%              SP-42, pp1038-1051, 1994

% ------- dttmtxno.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 the forward DCT type II even
if nargin < 2; type = 'dct2e'; end
if nargin < 3; frwd = 1;       end
if nargin < 4; 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/2;
l      = ones(N,1);
l(end) = 1/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 = 2    *cos( m     .* n     *pi /N) *k;
    case {'dctiie', 'dct2e'}
        T = 2    *cos( m     .*(n+.5) *pi /N);
    case {'dctiiie','dct3e'}
        T = 2    *cos((m+.5) .* n     *pi /N) *k;
    case {'dctive', 'dct4e'}
        T = 2    *cos((m+.5) .*(n+.5) *pi /N);
    % ----------------------------------------------------
    % DCT I-IV odd
    case {'dctio',  'dct1o'}
        T = 2    *cos( m     .* n     *2 *pi /(2*N-1)) *k;
    case {'dctiio', 'dct2o'}
        T = 2    *cos( m     .*(n+.5) *2 *pi /(2*N-1)) *l;
    case {'dctiiio','dct3o'}
        T = 2    *cos((m+.5) .* n     *2 *pi /(2*N-1)) *k;
    case {'dctivo', 'dct4o'}
        N = N+1;
        T = 2    *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 = 2    *sin( m     .* n     *pi /N);
    case {'dstiie', 'dst2e'}
        m = m+1;
        T = 2    *sin( m     .*(n+.5) *pi /N);
    case {'dstiiie','dst3e'}
        n = n+1; k(end) = k(1); k(1) = 1;
        T = 2    *sin((m+.5) .* n     *pi /N) *k;
    case {'dstive', 'dst4e'}
        T = 2    *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    *sin( m     .* n     *2 *pi /(2*N-1));
    case {'dstiio', 'dst2o'}
        m = m+1; N = N+1;
        T = 2    *sin( m     .*(n+.5) *2 *pi /(2*N-1));
    case {'dstiiio','dst3o'}
        n = n+1; N = N+1;
        T = 2    *sin((m+.5) .* n     *2 *pi /(2*N-1));
    case {'dstivo', 'dst4o'}
        T = 2    *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

% Calculate the scale factor for inverse
if ~frwd
    switch lower(type(end))
        case 'e'
            N = 2*N;
        case 'o'
            N = 2*N-1;
    end
    % Apply the scale factor
    T = T/N;
end