matlab - 使用一维 DCT 的二维离散余弦变换

Question

我正在尝试通过使用对图像实现二维离散余弦变换一维 DCT 运算。如果我将它与 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

Answer 1

您是正确的，因为 2D DCT 是可分离的。您只需先将 1D DCT 应用于每一行，然后获取中间结果并将其应用于列。但是，您有两个基本错误。让我们来看看它们。

错误 #1 - DCT 的大小不正确

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

l=size(signal,1);

因为您对每一行应用 DCT，然后对每一列应用 DCT，所以只有当您将 DCT 应用到列 时，以上内容才有效。 size(signal,1) 如果输入是一列，肯定会给出输入向量的长度。但是，如果您的输入是行，则size(signal,1) 的输出将为1。因此，您应该将 size(signal,1) 替换为 numel，这样您就可以确保获得元素的总数 - 无论输入是一行还是专栏。

此外，如果您想让代码兼容以在 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

---

如您所见，这两个结果是匹配的。我觉得不错!

---

为了确保我们的一致性，这是您的两个函数的完整代码 list ，以及我所做的修改:

% 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