离散余弦变换一维Matlab

Discrete Cosine Transformation 1D Matlab

我正在尝试使用矩阵向量乘法在 MATLAB 中实现 DCT。具体来说,我想创建系数的 DCT 矩阵,然后用它与一维信号相乘来计算一维 DCT。

这是我生成 DCT 矩阵的代码:

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

在此之后,我尝试使用 x = [1 2 3 4 5 6 7 8]:

的测试信号执行 DCT
x=[1 2 3 4 5 6 7 8];
y=dct1(8)*x';

它给出的答案是:

12.7279
18.0000
18.0000
18.0000
18.0000
18.0000
18.0000
18.0000

然而,正确答案是:

12.7279
-6.4423
-0.0000
-0.6735
      0
-0.2009
-0.0000
-0.0507

您的代码中的错误非常轻微但很重要。您计算系数的行:

D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);

看看这行的最后一部分:

D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/2*n);
%//                              ^^^

Because multiplication and division are equal in precedence,这和做完全一样:

D(u+1,x+1)=au*cos((((2*x+1)*u*pi)/2)*n);

因此,您没有除以 2n。你除以 2 然后 乘以 n 这是不正确的。您只需要用方括号将 2*n 操作括起来:

D(u+1,x+1)=au*cos(((2*x+1)*u*pi)/(2*n)); 

一旦你这样做了,我们就得到了正确的 DCT 矩阵。顺便说一句,检查您是否有正确答案的一种方法是使用 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,然后计算元素乘积。

这样的东西会起作用:

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

不是使用 if 语句来确定每行需要应用的权重,我们可以只使用 sqrt(2/n) 然后除以第一行的 sqrt(2) 这样你而是除以 sqrt(1/n)。此代码应产生与更正代码相同的结果1.


无论如何,一旦我做了这些更正,我就比较了你的代码给出的答案和 dctmtx 给出的答案,它是正确的:

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

一旦我们得到校正后的 DCT 矩阵,我们可以通过与您使用的 1:8 的测试向量进行矩阵乘法来检查实际的 DCT 计算:

>> dct1d(8)*((1:8).')

ans =

   12.7279
   -6.4423
   -0.0000
   -0.6735
         0
   -0.2009
   -0.0000
   -0.0507

1.这段代码实际上是在 dctmtx 中所做的,一旦你去掉所有的错误检查和输入一致性检查。