% 2d discrete cosine transform matrix: constructor
%
% A = matDCT2(n,N)
%
% Function does not rely on Matlab's Signal Processing Toolbox.
%
%  n  [1,1]: side length of the quadratic image to transform
%  N  [1,1]: number of components to keep [1,..,n^2], default N = n^2
%  complex [1,2]: complex number treatment A:IN^n->OUT^m,          default [1,1]
%                 first number specifies output field OUT^m and the second the
%                 input field IN^m of the matrix A
%                 1 -  real(z) + 1i*imag(z)    IC:   native complex numbers
%                 2 - [real(z); imag(z)]       IR^2: stacked real numbers
%                 3    real(z)                 IR:   real part only
%
% (c) by Hannes Nickisch, MPI for Biological Cybernetics, 2010 August 23

function A = matDCT2(n,N,complex)

  if nargin < 1, n = 1;   end, n = fix(n); if n<1, n = 1; end
  if nargin < 2, N = n*n; end, N = fix(N); if N<1, N = 1; end
  if N>n^2, error('N cannot be greater than n^2'), end

  k=1; id = zeros(n,n);                  % natural index order in the DCT domain
  for i=1:n
    for j =1:i,      id(i, j) = k; k=k+1; end
    for i2=i-1:-1:1, id(i2,j) = k; k=k+1; end
  end
  [ord,ord] = sort(id(:),'ascend');         % prior order of cosine coefficients
  A.coeff = ord(1:N);
  sz = [N,n^2]; % matrix size
  A.imsz = [n,n]; % image size
  if nargin<3, complex = 1; end
  if numel(complex)==1, complex = complex*[1,1]; end
  A = class(A,mfilename,mat(sz,complex));
