
function C = dctallmtx(N,type)
% DCTALLMTX Discrete cosine transform matrix
%   C = dctallmtx(N,type)
%
%   This generates all the DCT transformation matrices I-IV.
%   Type II is the most widely used.
%
%   N:    transform size (NxN)
%   type: dct type, (1:4)
%   C:    the transformation matrix
%
%   Notes:  The implementation generates unitary DCT matrices.
%   Here are some interesting properties
%   The inverse of types I and IV is itself (involutary)
%   The inverse of type II is type III and vice versa
%   The inverse of type II is its transpose
%   The inverse of type III is its transpose
%
%   Example:
%   Transforming an framed input signal x (frame per column)
%   can be done as:
%   [N,M] = size(x);
%   C = dctmtx(N,1); % Type 1
%   X = C*x;
%
%   Reference:  "Discrete Cosine Transform" K.R. rao, P. Yip
%               ISBN 0-12-580203-X,  pp 10-16

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

% The default type is II
if nargin < 2; type = 2; end

% K is the size of the transform
K = N;

% For Type I change the N (this is different than pp.11)
if type == 1; N = N-1; end

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

% ... and the scaling factors
p       = sqrt(2/N);
k       = ones(K,1);
k(1)    = 1/sqrt(2);
if type == 1; k(end) = k(1); end

% Make k a diagonal matrix 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);

% Make the transformation matrix
switch type
    case 1
        C = p *k *cos( m     .* n     *pi /N) *k;
    case 2
        C = p *k *cos( m     .*(n+.5) *pi /N);
    case 3
        C = p    *cos((m+.5) .* n     *pi /N) *k;
    case 4
        C = p    *cos((m+.5) .*(n+.5) *pi /N);
    otherwise
        error('Unknown DCT type, try 1:4')
end
