idctt(b,nn,omega,FFTW)

function a = idctt(b,nn,omega,FFTW)
%IDCTT Inverse discrete cosine transform
%   Specifically, a modified Type III DCT --
%   it computes an orthogonal DCT, 
%   i.e. IDCT = (DCT)^* = (DCT)^-1, where the 
%   foward DCT is modified Type II
%   Unless give the fftw argument (see below), this gives the same
%   result as calling matlab's IDCT
%
%   modified from IDCT (builtin matlab function) by Stephen Becker
%
%   X = IDCTT(Y) inverts the DCT transform, returning the
%   original vector if Y was obtained using Y = DCT(X).
%
%   X = IDCTT(Y,N) pads or truncates the vector Y to length N 
%   before transforming.
%
%   Y = IDCTT(X,N,OMEGA) uses Y = IDCTT(XX), where XX is a zero
%   vector (or matrix) of size N, and XX(OMEGA) = X
%   This is the inverse of DCTT(X,N,OMEGA), which is a sampled DCT.
%   You MUST supply the argument N for this case
%
%   Y = IDCTT(X,N,OMEGA,1) will return a Type III DCT using
%   the normalization coefficients from the fftw website
%   (and 2x the coefficients from the wikipedia definitions)
%   This version is NOT orthogonal.  It computes an
%   inverse to the forward DCT of fftw type II up to scaling, i.e.
%       DCT(IDCT) = IDCT(DCT) = 2n * I
%   where DCT and IDCT are DCTT and IDCTT, resp., with the fftw flag.
%   This is what fftw will compute (if you specify type III)
%   DCT and IDCT are no longer adjoints, i.e. DCT no longer othogonal.
%
%   Y = IDCTT(X,N,OMEGA,2) will return a Type III DCT
%   that is the ADJOINT of the forward DCT of fftw type II, i.e.
%       DCT = (IDCT)^*
%   Note that 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 Y is a matrix, the IDCT operation is applied to
%   each column.
%
%   See also DCTT.

%   
%   This version is modified by Stephen Becker, 1/28/08
%   [e-mail]
%   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 IDCT
%
%   Based off matlab's "idct" 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

if nargin == 0,
	error('Not enough input arguments.');
end

if isempty(b),
   a = [];
   return
end

% If input is a vector, make it a column:
do_trans = (size(b,1) == 1);
if do_trans, b = b(:); end
   
[n,m] = size(b);
if nargin==1 || isempty(nn)
  nn = size(b,1);
end

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

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

% Pad or truncate b if necessary
if n<nn,
  bb = zeros(nn,m);
  if SUBSAMPLE
      bb(omega,:) = b;
  else
      bb(1:n,:) = b;
  end
elseif n==nn
  bb = b;
else
  bb = b(1:nn,:);
end
n = nn;

if nargin > 3 && FFTW
% If requested, convert to fftw conventions
%   FFTW = 1 computes the Inverse DCT, up to a constant of 2*n
%   FFTW = 2 computes the Adjoint DCT
% For all except the DC component (k=0),
%   multiply by sqrt(n/2)
% For the DC component, multiply by just sqrt(n)
    if FFTW == 2   sqrtn = 2*sqrt(n);
    else           sqrtn =   sqrt(n);
    end
    sqrtn2 = sqrt(2*n);
    if m == 1   % break into subcases for speed
        bb(2:end) = sqrtn2 * bb(2:end);
        bb(1) = sqrtn * bb(1);
    else
        bb(2:end,:) = sqrtn2 * bb(2:end,:);
        bb(1,:) = sqrtn * bb(1,:);
    end
end

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

% Compute weights with persistent variables
if RECOMPUTE
    c = j*pi/(2*n);
    ww = sqrt(2*n) * exp( c*(0:n-1) ).';
    if odd
       ww(1) = ww(1) * sqrt(2);
    else
       ww(1) = ww(1)/sqrt(2);
    end
end

if odd
  % Form intermediate even-symmetric matrix.
  if m == 1
      W = ww; % matlab is smart enough NOT to perform a copy here
  else
      W = ww(:,ones(1,m));   % MAKE THIS INTO A REPMAT CALL?
      % Note to user: if you repeatedly call this function
      % with a *matrix* argument "b", then you might want to make
      % "W", not just "w", a persistent variable, if you can afford
      % the memory.
  end
  
  % most of the code below here is identical to the matlab version
  yy = zeros(2*n,m);
  yy(1:n,:) = W.*bb;
  yy(n+2:2*n,:) = -j*W(2:n,:).*flipud(bb(2:n,:));
  
  y = ifft(yy);

  % Extract inverse DCT
  a = y(1:n,:);

else
  if m == 1
      yy = ww.*bb;
  else
      W = ww(:,ones(1,m));
      yy = W.*bb;
  end

  % Compute x tilde using equation (5.93) in Jain
  y = ifft(yy);
  
  % Re-order elements of each column according to equations (5.93) and
  % (5.94) in Jain
  a = zeros(n,m);
  a(1:2:n,:) = y(1:n/2,:);
  a(2:2:n,:) = y(n:-1:n/2+1,:);
end

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