1
$\begingroup$

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;
Pixel=imread('small.png');
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);
$\endgroup$
  • 1
    $\begingroup$ Have you compared the results to a known good implementation? Matlab certainly has DCT functions. $\endgroup$ – Hilmar Sep 30 at 12:30
  • 1
    $\begingroup$ 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. $\endgroup$ – Florian Sep 30 at 13:36
  • $\begingroup$ 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. $\endgroup$ – Suvi Sep 30 at 15:44
  • $\begingroup$ @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? $\endgroup$ – Suvi Sep 30 at 16:27
  • $\begingroup$ 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 $\endgroup$ – Juha P Sep 30 at 17:13
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Suvi Oct 2 at 16:38
  • $\begingroup$ 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). $\endgroup$ – Florian Oct 2 at 21:06
  • $\begingroup$ But if I exchange n' and q to n and q', then C0 becomes a 1-by-1-matrix. $\endgroup$ – Suvi Oct 5 at 10:47
  • $\begingroup$ Well yeah, you need to exchange the order too. I edited the reply. $\endgroup$ – Florian Oct 5 at 12:10
  • $\begingroup$ 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. $\endgroup$ – Suvi Oct 6 at 10:11
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you, @Juha P. I got it installed. $\endgroup$ – Suvi Oct 2 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.