Code covered by the BSD License

# idctt

### Stephen Becker (view profile)

28 Feb 2008 (Updated )

type II Inverse DCT (i.e. type III DCT)

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.
%

%
%   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 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```