I'm trying to implement DCT in MATLAB by using matrix-vector multiplication. Specifically, I'd like to create the DCT matrix of coefficients then use this to multiply with 1D signal to compute the 1D DCT.
Here's my code to produce the DCT matrix:
function D=dct1d(n)
for u=0:n-1
if u==0
au=sqrt(1/n);
else
au=sqrt(2/n);
end
for x=0:n-1
D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);
end
end
After this I tried performing the DCT with a test signal of x = [1 2 3 4 5 6 7 8]:
x=[1 2 3 4 5 6 7 8];
y=dct1(8)*x';
The answer it gives is:
12.7279
18.0000
18.0000
18.0000
18.0000
18.0000
18.0000
18.0000
However, the correct answer is:
12.7279
-6.4423
-0.0000
-0.6735
0
-0.2009
-0.0000
-0.0507
The error in your code is very slight but fundamental. The line where you're computing the coefficients:
D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);
Take a look at the very last part of the line:
D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);
%// ^^^
Because multiplication and division are equal in precedence, this is exactly the same as doing:
D(u+1,x+1)=au*cos((((2*x+1)*u*pi)/2)*n);
As such, you aren't dividing by 2n. You are dividing by 2 then multiplying by n which is incorrect. You simply have to surround the 2*n operation with brackets:
D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/(2*n));
Once you do this, we get the correct DCT matrix. BTW, one way to check if you have the right answer is to use the dctmtx function, which calculates the N x N DCT coefficient matrix that you're seeking. However, this function is part of the Signal Processing Toolbox, so if you don't have that then you unfortunately can't use this function, but if I can suggest an alternative answer instead of using for loops, I would build a 2D grid of coordinates with meshgrid, then compute element-wise products.
Something like this would work instead:
function D = dct1d(n)
[x,u] = meshgrid(0:n-1);
D = sqrt(2/n)*cos(((2*x+1).*u*pi)/(2*n));
D(1,:) = D(1,:) / sqrt(2);
end
Instead of using an if statement to determine what weights we need to apply per row, we can just use sqrt(2/n) then divide by sqrt(2) for the first row so that you're dividing by sqrt(1/n) instead. This code should produce the same results as your corrected code1.
In any case, once I made these corrections, I compared both of the answers between what your code gives and what dctmtx gives and it's right:
>> dct1d(8)
ans =
0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536
0.4904 0.4157 0.2778 0.0975 -0.0975 -0.2778 -0.4157 -0.4904
0.4619 0.1913 -0.1913 -0.4619 -0.4619 -0.1913 0.1913 0.4619
0.4157 -0.0975 -0.4904 -0.2778 0.2778 0.4904 0.0975 -0.4157
0.3536 -0.3536 -0.3536 0.3536 0.3536 -0.3536 -0.3536 0.3536
0.2778 -0.4904 0.0975 0.4157 -0.4157 -0.0975 0.4904 -0.2778
0.1913 -0.4619 0.4619 -0.1913 -0.1913 0.4619 -0.4619 0.1913
0.0975 -0.2778 0.4157 -0.4904 0.4904 -0.4157 0.2778 -0.0975
>> dctmtx(8)
ans =
0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536
0.4904 0.4157 0.2778 0.0975 -0.0975 -0.2778 -0.4157 -0.4904
0.4619 0.1913 -0.1913 -0.4619 -0.4619 -0.1913 0.1913 0.4619
0.4157 -0.0975 -0.4904 -0.2778 0.2778 0.4904 0.0975 -0.4157
0.3536 -0.3536 -0.3536 0.3536 0.3536 -0.3536 -0.3536 0.3536
0.2778 -0.4904 0.0975 0.4157 -0.4157 -0.0975 0.4904 -0.2778
0.1913 -0.4619 0.4619 -0.1913 -0.1913 0.4619 -0.4619 0.1913
0.0975 -0.2778 0.4157 -0.4904 0.4904 -0.4157 0.2778 -0.0975
Once we get the corrected DCT matrix, we can check the actual DCT calculations by performing the matrix multiplication with the test vector of 1:8 that you used:
>> dct1d(8)*((1:8).')
ans =
12.7279
-6.4423
-0.0000
-0.6735
0
-0.2009
-0.0000
-0.0507
1. This code is actually what is being done in dctmtx under the hood, once you strip away all of the error checking and input consistency checks.