# DCT with Matlab

I tried to create a Matlab/Octave code to implement this DCT formula.

Is the code below correct? As a help I used the code from this site and modified it a bit

clc
clear
close all
N=8;
figure
imshow(Pixel)
Cosines=zeros(N,N);
for x1=1:1:N
for i1=1:1:N
Cosines(x1,i1)=cos( ( 2*pi*x1*(i1+0.5) )/ (2*N) );
endfor
endfor
Coeffi=zeros(N,1)+sqrt(2);
Coeffi(1)=1;
Coefficient=zeros(N,N);
for x=1:1:N
for y=1:1:N
Coefficient(x,y)=Coeffi(x)*Coeffi(y);
endfor
endfor
for i=1:1:N
for j=1:1:N
temp = 0.0;
for x=1:1:N
for y=1:1:N
temp = temp+(Cosines(x,i)*Cosines(y,j)*Pixel(x,y));
endfor
endfor
temp = (1/N) * Coefficient(i,j);
DCT(i,j) = temp;
endfor
endfor
figure
imshow(DCT);

• Have you compared the results to a known good implementation? Matlab certainly has DCT functions. – Hilmar Sep 30 at 12:30
• I was wondering too: any reasons you haven't used the dct2 function? Otherwise: you can get rid of all the loops in your script by using matrices-vector products. This is not only easier to read but also much faster. – Florian Sep 30 at 13:36
• The reasons why I don't use the DCT function are: - I have only Octave, not Matlab (Octave has no DCT functions - This is for learning purpose. I need to understand the implementation, not just use it. – Suvi Sep 30 at 15:44
• @Florian How would you write it using matrices-vector products? I tried this [dropbox.com/s/r4ksa1wj9boo3ln/dctcode.m?dl=0], what do you think about it? – Suvi Sep 30 at 16:27
• Moi, are you using latest version of Signal package octave.sourceforge.io/signal/index.html ? Signal package lists DCT and DCT2 : octave.sourceforge.io/signal/function/dct.html , octave.sourceforge.io/signal/function/dct2.html – Juha P Sep 30 at 17:13

As you are asking how to do this without loops, here is what I would do:

N = 8;
n = (0:N-1);
q = (0:N-1);
C0 = cos(pi*(2*n'+1)*q/2/N);                 % unscaled DCT matrix
C = C0 * sqrt(diag([1, 2*ones(1,N-1)])/N);   % column scaling

% let's try if that is consistent with the regular DCT
v = rand(N,1);
disp(norm(C*v - dct(v,'type',3))); % it's actually a "type 3", there are many DCTs

% now a 2-D DCT looks like this
Image = rand(N,N);
Image_DCT = C*Image*C';


The code uses the fact that the product of a column vector (n') and a row vector (q) gives a rank one matrix of all pairwise products. This should work even in old versions of Matlab or Octave that have no broadcasting functionality yet. Give it a try.

Not sure this is exactly the DCT you need as there are several versions of it. But you should be able to get them with some minor tweaks.

*edit: here is the version for a type-2 DCT

N = 8;
n = (0:N-1);
q = (0:N-1);
C0 = cos(pi*q'*(2*n+1)/2/N);                 % unscaled DCT matrix
C = sqrt(diag([1, 2*ones(1,N-1)])/N) * C0;   % row scaling

% let's try if that is consistent with the regular DCT
v = rand(N,1);
disp(norm(C*v - dct(v,'type',2))); % it's actually a "type 2", there are many DCTs

% now a 2-D DCT looks like this
Image = rand(N,N);
Image_DCT = C*Image*C';


*edit2: As you asked about the mathematics, if we consider the 2-D DCT as a matrix-matrix product in the form $$D = C\cdot I\cdot C^T$$ then its $$(p,q)$$ element is given by $$D_{p,q} = \sum_m \sum_n C_{p,m} I_{m,n} C_{q,n}.$$ Since each of the $$C_{i,j}$$ represents one of the cosine functions, this should explain how this computation aligns with the scalar representation you saw in the video.

• Thank you, @Florian. I compared your code with the dctmtx() function. Your code produces a transformation matrix (C) which is the transpose of what dctmtx() produces. – Suvi Oct 2 at 16:38
• Yeah, if you need the transpose (dct calls it 'Type 2') you can just transpose C. Alternatively, exchange n' and q to n and q', respectively (and change column to row scaling by pre- instead of postmultiplying with the sqrt(diag(...)) matrix). – Florian Oct 2 at 21:06
• But if I exchange n' and q to n and q', then C0 becomes a 1-by-1-matrix. – Suvi Oct 5 at 10:47
• Well yeah, you need to exchange the order too. I edited the reply. – Florian Oct 5 at 12:10
• Ok thanks. Can somebody explain me why for example in this video (youtube.com/watch?v=DS8N8cFVd-E) there is cos() times cos() although in Florian's code there is CImageC' so cos is not multiplied with cos but with the pixels of the image. So it's not C*C' but there's A between the cosines. – Suvi Oct 6 at 10:11

You can test your code also in Octave (v. 4.0 and above) by first installing two packages for Octave; Signal - https://octave.sourceforge.io/signal/index.html Control - https://octave.sourceforge.io/control/index.html

DCT and DCT2 functions are parts of Signal.

Package installation instructions can be found from here.

• Thank you, @Juha P. I got it installed. – Suvi Oct 2 at 16:39