1D DCT matlab code

Question

I was writing MATLAB code to compute 1D DCT of sample y. On computing DCT for y=[0,1,2], code generates coefficient [3.0000 -2.2304 0 -0.1585] which was consistent to coefficient generated by Matlab default dct([0,1,2]).

When sample y=[1,2,3] was chosen the same code didn't work well.

Looking forward to modification and clarification in my custom DCT code.

Eq DCT:

$$C(u) = \alpha(u)\sum_{x=0}^{N-1}f(x)\cos\left[\dfrac{\pi(2x+1)u}{2N} \right] \qquad\text{for}\qquad u=0,1,2,\ldots,N-1$$

Custom code

y=[0,1,2];
[M,N]=size(y);
sum=0;
dct1d=zeros(1,N);
u=[0:N-1] ; 
for j=1:N
    for i=1:N
        sum=sum+y(i).*(cos((pi.*(2.*y(i)+1).*u(j))/(2*N)));
    end
    if j==1   
            K=sqrt(1/N);
    else
            K=sqrt(2/N);
    end

   dct1d(j)=K.*sum;
   sum=0;
end
dct1d

Answer 1

You have mistyped the formula, replace this line

sum = sum + y(i).*(cos((pi.*(2.*y(i)+1).*u(j))/(2*N)));

with the one below, and it works fine.

sum = sum + y(i).*(cos((pi.*(2.*u(i)+1).*u(j))/(2*N)));

Answer 2

For clarity, I would write this DCT as:

$$F(u) = \alpha(u)\sum_{i=0}^{N-1}f(i)\cos\left(\frac{\pi u}{2N}(2i+1)\right)$$ We note that, with this 1-indexing of Matlab:

$$y[i+1] = f(i)\,.$$

Then I would modify the inner limit (from $0$ to $N-1$ instead of $1$ to $N$):

for i = 0:N-1
sum = sum + y(i+1)*(cos((pi*(2*i+1)*u(j))/(2*N)));
end

and you can remove the dotted operator as well. I do prefer keeping the discrete time/space index $i$ better separated from the discrete frequency index $u$, to preserve the symmetry/duality of the Discrete Cosine Transform, and that of its inverse as well. This could be beneficial for extensions or generalizations of similar trigonometric transforms.