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