dct.m, anyone?

[Clay M. Thompson]

In article <1993Feb8.1...@bert.eecs.uic.edu> Steven L. Eddins,
edd...@bert.eecs.uic.edu writes:
>Does anyone have an m-file to compute the DCT? Yes, I know there are
>several definitions; any would be acceptable.

Here are files for dct and idct that we have at The MathWorks. Enjoy!
By the way, these functions work with MATLAB 4.0 (Hint: size(a,1) is
not compatible with MATLAB 3.5).

- Clay (cl...@mathworks.com)

---------------------------------
function b=dct(a,n)
%DCT Discrete cosine transform.
%
% Y = DCT(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 = DCT(X,N) pads the vector X to length N before
% transforming.
%
% If X is a matrix, the DCT operation is applied to each
% column. This transform can be inverted using IDCT.
%
% See also: FFT,IFFT, and IDCT.

% Clay M. Thompson 2-12-93
% Copyright (c) 1993 by The MathWorks, Inc.

if nargin==1,
if min(size(a))==1, a = a(:); end
n = size(a,1);
end
m = size(a,2);

% Pad a if necessary
if size(a,1)<n,
aa = zeros(n,m);
aa(1:size(a,1),:) = a;
else
aa = a;
end

% Form intermediate even-symmetric matrix.
y = zeros(2*n,m);
y(1:n,:) = aa;
y(n+1:n+n,:) = flipud(aa);

% Perform FFT
yy = fft(y);

% Compute DCT coefficients
W = exp(((0:n-1)'* (ones(1,m)*(-sqrt(-1)*pi/2/n))));
b = W.*yy(1:n,:);

if ~any(any(imag(a))), b = real(b); end

----------------------
function a = idct(b,m,n)
%IDCT Inverse discrete cosine transform.
%
% X = IDCT(Y) inverts the DCT transform, returning the
% original vector if Y was obtained using Y = DCT(X).
%
% X = IDCT(Y,N) pads the vector Y to length N before
% transforming.
%
% If Y is a matrix, the IDCT operation is applied to
% each column.
%
% See also: FFT,IFFT,DCT.

% Clay M. Thompson 2-12-93
% Copyright (c) 1993 by The MathWorks, Inc.

if nargin==1,
if min(size(b))==1, b = b(:); end
n = size(b,1);
end
m = size(b,2);

% Pad b if necessary
if size(b,1)<n,
bb = zeros(n,m);
bb(1:size(b,1),:) = b;
else
bb = b;;
end

% Form intermediate even-symmetric matrix.
W = exp(([0:n-1]'* (ones(1,m)*(sqrt(-1)*pi/2/n))));
yy = zeros(2*n,m);
yy(1:n,:) = W.*bb;
yy(n+2:n+n,:) = -sqrt(-1)*W(2:n,:).*flipud(bb(2:n,:));

y = ifft(yy);

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

if ~any(any(imag(b))), a = real(a); end

[Saad John Bedros]

I extracted this code from a larger file
it includes the dct fuction which I am hoping to mex it
one of these days

%load boat128 ; im = boat128 ; clear boat16 ;
%end
cl = 2 ;
class = 2 .^ cl
%read the image
[imsizex ,imsizey] = size(im) ;

imim = zeros(imsizex ,imsizey) ;
N = 8 ;
NN = N*N ;
vectnum = imsizex * imsizey / NN ;

dctim = zeros(vectnum,NN) ;
dctf = zeros(N,N) ;
fkrnl = dctkernel(N) ;

t0 = clock ;
k = 1 ;
for i=1:N:imsizex,
for j=1:N:imsizey,
dctim(k,:) = rect2zigzag(fkrnl* (im(i:(i+N-1),j:(j+N-1))-128) *fkrnl') ;
k = k+1 ;
end
end

% DCT KERNEL
function [krnl] = dctkernel (NN)

% Calculating DCT constant of 1-D.
AU = ones(1,NN) /NN ;
AU (1) = AU(1)/sqrt(2.);
AU = 2 * AU ;

NN2 = 1 / (2 * NN);
% Calculating DCT kernel of 1-D.
u = 0:NN-1 ;
x = u ;
upi = (u * pi) * NN2 ;
cs = cos (upi' * (2*x+1) ) ;
krnl = (AU' * ones(1,NN)) .* cs;

Saad Bedros
Elect Eng Dept
U of Minnesota