离散余弦变换1D Matlab

Question

离散余弦变换1D Matlab

matlabsignal-processingimage-compressiondct

3

我试图在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
-6.4423
-0.0000
-0.6735
      0
-0.2009
-0.0000
-0.0507

- Mark Ben

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- rayryeng · Accepted Answer

你的代码中出现了一个非常微小但基本的错误。在计算系数的那一行：

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);
%//                              ^^^

由于乘法和除法具有相同的优先级，因此这与执行以下操作完全相同：

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

不需要使用 if 语句来确定每行需要应用的权重，我们可以只使用 sqrt(2/n)，然后在第一行除以 sqrt(2)，这样你就是在除以 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 在幕后所执行的操作。}