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

score 4 · Accepted Answer · answered Jun 30 '20 at 14:32

4

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)));

answered Jun 30 '20 at 14:32

Fat32

28,152
3
24
50

I voted for you, yet I am wondering whether it is a good idea to multiple a $u(i)$ by some $u(j)$. Same result here, but we lose the idea that inside the cosine, there is a time index and a frequency index (duality). Combining them cause issues with generalized versions (ie with different offsets in both index lists) – Laurent Duval Jun 30 '20 at 15:19
1

@LaurentDuval His style is his convention. However I can see that the integers i and j refers to time and frequency indices. He uses an array u for index remapping which I found as a clever thing... – Fat32 Jul 01 '20 at 19:46
The problem was to have $y$ both inside and outside the $\cos$. The philosophical question is: a unitless remapping for variables, or two distinct remappings with inverse units and same values – Laurent Duval Jul 01 '20 at 21:14

Laurent Duval · Answer 2 · 2020-06-30T18:54:07.840

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.

1D DCT matlab code

2 Answers2

Linked