使用1D DCT的2D离散余弦变换

Question

我正在尝试使用1D DCT操作将2D离散余弦变换应用于图像。如果我将其与dct2 MATLAB函数进行比较，我的输出结果是不正确的。我不明白我的代码出了什么问题以及在哪里出了问题。

如果有人能指出错误或提供任何其他建议，那将非常有帮助。

这是我用MATLAB编写的代码：

% main function
signal=rand(100);
signal_dct=myDCT(signal);
figure; imshow((signal_dct));

% function to calculate 2D DCT of an image
function res=myDCT(signal)

    signal=double(signal);

    l=size(signal,1);

    res=zeros(l);  %initialize the final result matrix

    for k=1:l     %calculate 1D DCT of each row of image
        res(k,:)=mdct(signal(k,:));  
    end

    for k=1:l   %calculate 1D DCT of each column of image
        res(:,k)=mdct(res(:,k));
    end
end

%% function to calculate 1D DFT of a 1D signal
function res=mdct(signal)

    l=size(signal,1);

    for i=1:l

        if i==1   %for signal index of 1, alpha is 1/sqrt(l)
            alpha=sqrt(1/l);
        else   %for signal index of greater than 1
            alpha=sqrt(2/l);
        end

        j=[1:l];
        % summation calculates single entry of res by applying the  
        % formula of DCT on the signal
        summation=sum(sum(signal(j)*cos((pi*(2*(j-1)+1)*(i-1))/(2*l))));
        res(i)=alpha*summation;        
    end
end

- bumpk_sync

---

看起来你正在使用以弧度为单位的cosd函数。请改用cos函数。 - Luis Mendo

谢谢你指出那个错误。我已经修正了它。但是结果仍然与我从" dct2 "Matlab内置函数得到的完全不同。 - bumpk_sync

我真的不知道。无论如何，如果您发布一个可重现的示例（可以运行的代码，并附有说明问题的实际数据），您更有可能得到答案。 - Luis Mendo

我已经包含了一些主要的代码行来查找随机二维数字数组的DCT。 - bumpk_sync

@rayryeng 做得好，回答很棒！ - Luis Mendo

Answer 1

您说得对，2D DCT是可分离的。您只需先将1D DCT应用于每一行，然后将中间结果应用于列即可。但是，您有两个基本错误。 让我们逐个解决。

错误1 - DCT的大小不正确

具体来说，请查看您的mdct函数中的此语句:

l=size(signal,1);

由于您正在对每行，然后每列应用DCT，因此如果您正在将DCT应用于列，则上面的方法只适用于该情况。如果输入是一列，则size(signal,1)肯定会为您提供输入向量的长度。但是，如果您的输入是一行，则size(signal,1)的输出将为1。因此，您应该使用numel替换size(signal,1)，以确保您将获得元素的总数 - 无论输入是行还是列。

另外，如果您想使代码兼容在DCT循环中执行求和操作，您应该始终确保输入是行向量。因此，请改用以下方式进行操作：

l = numel(signal);
signal = signal(:).';

第一行确定我们的输入信号有多少元素，第二行确保我们有一个行向量。这是通过(:)将元素展开为列向量，然后在之后添加.'来确保我们对结果进行转置以获得行向量。

错误＃2-求和语句不正确

接下来，您需要在求和中进行逐元素乘法，以获得您要查找的内容。您还不需要那个额外的sum调用。 这是多余的。因此，请将您的求和语句修改为：

summation=sum(signal.*cos((pi*(2*(j-1)+1).*(i-1))/(2*l)));

不需要执行signal(j)，因为j跨越了整个向量的长度，所以你只需要使用signal即可。

---

一旦我进行了这些更改，并且我在较小的矩阵上执行以确保我们获得相同的结果：

rng(123123);
signal=rand(7);
signal_dct=myDCT(signal);
signal_dct2 = dct2(signal);

最后一行代码调用dct2函数，以便我们可以比较您自定义函数的结果和dct2给我们的结果。

我们得到：

>> signal_dct

signal_dct =

    3.7455   -0.1854   -0.1552    0.3949    0.2182   -0.3707    0.2621
   -0.2747    0.1566   -0.0955    0.1415    0.3156   -0.0503    0.8581
   -0.2095    0.0233   -0.2769   -0.4341   -0.1639    0.3700   -0.2282
   -0.0282    0.0791    0.0517    0.4749   -0.0169   -0.4327    0.0427
   -0.4047   -0.4383    0.3415   -0.1120   -0.0229    0.0310    0.3767
   -0.6058   -0.0389   -0.3460    0.2732   -0.2395   -0.2961    0.1789
   -0.0648   -0.3173   -0.0584   -0.3461   -0.1866    0.0301    0.2710

>> signal_dct2

signal_dct2 =

    3.7455   -0.1854   -0.1552    0.3949    0.2182   -0.3707    0.2621
   -0.2747    0.1566   -0.0955    0.1415    0.3156   -0.0503    0.8581
   -0.2095    0.0233   -0.2769   -0.4341   -0.1639    0.3700   -0.2282
   -0.0282    0.0791    0.0517    0.4749   -0.0169   -0.4327    0.0427
   -0.4047   -0.4383    0.3415   -0.1120   -0.0229    0.0310    0.3767
   -0.6058   -0.0389   -0.3460    0.2732   -0.2395   -0.2961    0.1789
   -0.0648   -0.3173   -0.0584   -0.3461   -0.1866    0.0301    0.2710

---

正如您所看到的，这两个结果都相符。对我来说看起来很好！

---

为了确保我们是一致的，这是您的两个函数的完整代码清单，带有我所做的修改：

% function to calculate 2D DCT of an image
function res=myDCT(signal)

    signal=double(signal);

    l=size(signal,1);

    res = zeros(l);

    for k=1:l     %calculate 1D DCT of each row of image
        res(k,:)=mdct(signal(k,:));  
    end

    for k=1:l   %calculate 1D DCT of each column of image
        res(:,k)=mdct(res(:,k));
    end
end

%% function to calculate 1D DFT of a 1D signal
function res=mdct(signal)

    %// Change
    l = numel(signal);
    signal = signal(:).';

    for i=1:l

        if i==1   %for signal index of 1, alpha is 1/sqrt(l)
            alpha=sqrt(1/l);
        else   %for signal index of greater than 1
            alpha=sqrt(2/l);
        end

        j=[1:l];
        % summation calculates single entry of res by applying the  
        % formula of DCT on the signal
        %// Change
        summation=sum(signal.*cos((pi*(2*(j-1)+1).*(i-1))/(2*l)));
        res(i)=alpha*summation;        
    end
end

- rayryeng

---

非常感谢！我被卡住了，不知道出了什么问题。你指出的错误虽然很小，但在我的情况下可能会反复出现。下次我写Matlab代码时会记住它们的。我无法感谢你了。 :) - bumpk_sync

@bumpk_sync - 很荣幸:)谢谢您的接受。还请查看您的DFT帖子。我找出了问题并进行了纠正。 - rayryeng