离散余弦变换一维 Matlab

Question

我正在尝试通过使用矩阵向量乘法在 MATLAB 中实现 DCT。具体来说，我想创建系数的 DCT 矩阵，然后使用它与 1D 信号相乘以计算 1D 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

Answer 1

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

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 函数，该函数计算您正在寻找的N x N DCT 系数矩阵。但是，这个函数是信号处理工具箱的一部分，所以如果你没有那个，那么很遗憾你不能使用这个函数，但是如果我可以建议一个替代答案而不是使用 for 循环，我会建立一个带有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

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

---

无论如何，一旦我做出这些更正，我就会比较您的代码给出的答案和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 在幕后所做的，一旦你剥离了所有的错误检查和输入一致性检查。}