matlab - 2D Discrete Cosine Transform using 1D DCT -

i trying implement 2d discrete cosine transform image using 1d dct operations. output incorrect if compare dct2 matlab function. don't understand went wrong in code , it's happening.

if can point out mistake or other advice, helpful.

here code written in matlab

% main function signal=rand(100); signal_dct=mydct(signal); figure; imshow((signal_dct));  % function calculate 2d dct of image function res=mydct(signal)      signal=double(signal);      l=size(signal,1);      res=zeros(l);  %initialize final result matrix      k=1:l     %calculate 1d dct of each row of image         res(k,:)=mdct(signal(k,:));       end      k=1:l   %calculate 1d dct of each column of image         res(:,k)=mdct(res(:,k));     end end  %% function calculate 1d dft of 1d signal function res=mdct(signal)      l=size(signal,1);      i=1:l          if i==1   %for signal index of 1, alpha 1/sqrt(l)             alpha=sqrt(1/l);         else   %for signal index of greater 1             alpha=sqrt(2/l);         end          j=[1:l];         % summation calculates single entry of res applying           % formula of dct on signal         summation=sum(sum(signal(j)*cos((pi*(2*(j-1)+1)*(i-1))/(2*l))));         res(i)=alpha*summation;             end end 

you correct in 2d dct separable. apply 1d dct every row first, take intermediate result , apply columns. however, have 2 fundamental mistakes. let's go through them.

mistake #1 - size of dct isn't correct

specifically, @ statement in mdct function:

l=size(signal,1); 

because applying dct each row, each column, above work if applying dct columns. size(signal,1) give length of input vector if input column. however, if input row, output of size(signal,1) 1. therefore, should replace size(signal,1) numel sure going total number of elements - regardless if input row or column.

also, if want make code compatible summation in dct loop, should make sure input row vector regardless. such, instead:

l = numel(signal); signal = signal(:).'; 

the first line determines how many elements have our input signal, , second line ensures have row vector. done (:) unroll elements column vector, doing .' after ensure transpose result row vector.

mistake #2 - summation statement isn't correct

next, you're going have element-wise multiplications in summation you're looking for. don't need sum call there. it's superfluous. therefore, modify summation statement this:

summation=sum(signal.*cos((pi*(2*(j-1)+1).*(i-1))/(2*l))); 

there's no need signal(j) because j spans entire length of vector, , can signal.

---

once made these changes, , did on smaller size matrix ensure same results:

rng(123123); signal=rand(7); signal_dct=mydct(signal); signal_dct2 = dct2(signal); 

the last line of code calls dct2 can compare results custom function , dct2 gives us.

we get:

>> signal_dct  signal_dct =      3.7455   -0.1854   -0.1552    0.3949    0.2182   -0.3707    0.2621    -0.2747    0.1566   -0.0955    0.1415    0.3156   -0.0503    0.8581    -0.2095    0.0233   -0.2769   -0.4341   -0.1639    0.3700   -0.2282    -0.0282    0.0791    0.0517    0.4749   -0.0169   -0.4327    0.0427    -0.4047   -0.4383    0.3415   -0.1120   -0.0229    0.0310    0.3767    -0.6058   -0.0389   -0.3460    0.2732   -0.2395   -0.2961    0.1789    -0.0648   -0.3173   -0.0584   -0.3461   -0.1866    0.0301    0.2710  >> signal_dct2  signal_dct2 =      3.7455   -0.1854   -0.1552    0.3949    0.2182   -0.3707    0.2621    -0.2747    0.1566   -0.0955    0.1415    0.3156   -0.0503    0.8581    -0.2095    0.0233   -0.2769   -0.4341   -0.1639    0.3700   -0.2282    -0.0282    0.0791    0.0517    0.4749   -0.0169   -0.4327    0.0427    -0.4047   -0.4383    0.3415   -0.1120   -0.0229    0.0310    0.3767    -0.6058   -0.0389   -0.3460    0.2732   -0.2395   -0.2961    0.1789    -0.0648   -0.3173   -0.0584   -0.3461   -0.1866    0.0301    0.2710 

---

as can see, both results match up. looks me!

---

just sure consistent, full code listing both functions, modifications made:

% function calculate 2d dct of image function res=mydct(signal)      signal=double(signal);      l=size(signal,1);      res = zeros(l);      k=1:l     %calculate 1d dct of each row of image         res(k,:)=mdct(signal(k,:));       end      k=1:l   %calculate 1d dct of each column of image         res(:,k)=mdct(res(:,k));     end end  %% function calculate 1d dft of 1d signal function res=mdct(signal)      %// change     l = numel(signal);     signal = signal(:).';      i=1:l          if i==1   %for signal index of 1, alpha 1/sqrt(l)             alpha=sqrt(1/l);         else   %for signal index of greater 1             alpha=sqrt(2/l);         end          j=[1:l];         % summation calculates single entry of res applying           % formula of dct on signal         %// change         summation=sum(signal.*cos((pi*(2*(j-1)+1).*(i-1))/(2*l)));         res(i)=alpha*summation;             end end