function b=dctt(a,nn,omega,FFTW)
%DCTT  Discrete cosine transform of type II
%   Specifically, a modified Type II DCT --
%   it computes an orthogonal DCT, 
%   i.e. IDCT = (DCT)^* = (DCT)^-1, where the
%   IDCT is a DCT of modified Type III
%   Unless given the fftw argument (see below), this gives the same
%   result as calling matlab's DCT
%
%   modified from DCT (builtin matlab function) by Stephen Becker
%
%   Y = DCTT(X) returns the discrete cosine transform of X.
%   The vector Y is the same size as X and contains the
%   discrete cosine transform coefficients.
%
%   Y = DCTT(X,N) pads or truncates the vector X to length N 
%   before transforming.
%
%   Y = DCTT(X,N,OMEGA) returns Y(OMEGA), where OMEGA
%   is some subsample (i.e. so length(Y) = length(OMEGA) )
%   Use N = [] if you want to ignore this parameter
%   This is equivalent to sampling the rows of the DCT matrix.
%
%   Y = DCTT(X,N,OMEGA,true) will return a Type II DCT using
%   the normalization coefficients from the fftw website
%   (and 2x the coefficients from the wikipedia definitions).
%   This version is NOT orthogonal.  If the DCT and IDCT are
%   both called in this mode, then they are adjoints but not inverses,
%   and DCT(IDCT) is a diagonal matrix, and IDCT(DCT) is a full matrix.
%   Use N = [] and OMEGA = [] to ignore these parameters
%
%   If X is a matrix, the DCTT operation is applied to each
%   column.  This transform can be inverted using IDCTT.
%
%   See also IDCTT.
%
%   New behavior (6/19/08): if this function is called without any
%   arguments, the output is the number of function calls since
%   the last reset.  The calls are reset everytime this is queried.
%
%   Question: what are the FFTpack conventions?

%
%   This version is modified by Stephen Becker, 1/28/08
%   srbecker@caltech.edu
%   This version uses persistent variables to achieve
%   a more than 2x speed reduction on subsequent calls
%   On these subsequent calls, should be only about  
%   1.4 to 1.5x slower than a FFT,
%   instead of 4x slower (original matlab code)
%   On the first call, it's still faster than matlab's DCT
%
%   Based off matlab's "dct" function, which was written by:
%              C. Thompson, 2-12-93
%              S. Eddins, 10-26-94, revised

persistent ww   % this is independent of any subsampling
persistent N    
persistent nCalls
if isempty(nCalls) nCalls = 0; end

if nargin == 0,
    b = nCalls;
%     if nargout ==0 fprintf('Number of DCT calls is %d\n',b); end
    nCalls = 0;
    return
end
nCalls = nCalls + 1;
if isempty(a)
   b = [];
   return
end
% If input is a vector, make it a column:
do_trans = (size(a,1) == 1);
if do_trans, a = a(:); end

[n,m] = size(a);
if nargin<2 || isempty(nn)  % changing this slightly
  nn = size(a,1);
end
% m = size(a,2);

% decide whether to recompute "ww" (which depends only on n)
if isempty(N) || ( N ~= n)
    RECOMPUTE = true;
    N = n;
else
    RECOMPUTE = false;
end

% Pad or truncate input if necessary
if n<nn,
  aa = zeros(n,m);
  aa(1:n,:) = a;
elseif n==nn
  aa = a;
else
  aa = a(1:nn,:);    % truncate
end

odd =  ( rem(n,2)==1 || ~isreal(a) );

% Look at inputs and decide whether to subsample
if nargin > 2 && ~isempty(omega)
    SUBSAMPLE = true;
else
    SUBSAMPLE = false;
end

% Compute weights to multiply DFT coefficients
if RECOMPUTE
    c = -i*pi/(2*n);
    d = 1/sqrt(2*n);
    if ~odd, d = 2*d; end
    ww = d*(exp((0:n-1)*c)).';
    ww(1) = ww(1) / sqrt(2);
end
% ww = (exp(-i*(0:n-1)*pi/(2*n))/sqrt(2*n)).';  % SLOW
% ww(1) = ww(1) / sqrt(2);

if odd % odd case
  % Form intermediate even-symmetric matrix
  y = zeros(2*n,m);
  y(1:n,:) = aa;
  y(n+1:2*n,:) = flipud(aa);
  
  % Compute the FFT and keep the appropriate portion:
  yy = fft(y);  
  yy = yy(1:n,:);

else % even case
  % Re-order the elements of the columns of x
  if m == 1
      y = [ aa(1:2:n); aa(n:-2:2) ];
  else
      y = [ aa(1:2:n,:); aa(n:-2:2,:) ];
  end
  yy = fft(y);  
%   ww = 2*ww;  % Double the weights for even-length case  
end

% Multiply FFT by weights:
if m == 1
    if SUBSAMPLE
        b = ww(omega) .* yy(omega);
    else
        b = ww .* yy;
    end
else
    if SUBSAMPLE
        b = ww(omega,ones(1,m)) .* yy(omega,:);
    else
        b = ww(:,ones(1,m)) .* yy;
    end
end

if nargin > 3 && FFTW
% If requested, convert to fftw conventions
% For all except the DC component (k=0),
%   multiply by sqrt(n/2)
% For the DC component, multiply by just sqrt(n)
    sqrtn = 2*sqrt(n);
    sqrtn2 = sqrt(2*n);
    % Note that if we sample with omega, we might not have the DC component left
    if (SUBSAMPLE && (omega(1)==1)) || ~SUBSAMPLE
        if m == 1       % breaking into cases not necessary, but it's faster
            b(2:end) = sqrtn2*b(2:end);
            b(1)     = sqrtn *b(1);
        else
            b(2:end,:) = sqrtn2*b(2:end,:);
            b(1,:)     = sqrtn*b(1,:);
        end
    else
        b = sqrtn2*b;
    end
end

if isreal(a), b = real(b); end
if do_trans, b = b.'; end