使用1D DCT進行2D離散餘弦變換

Question

我正在嘗試使用1D DCT運算實現2D離散餘弦變換到圖像。如果我將它與MATLAB函數dct2進行比較，我的輸出是不正確的。我不明白我的代碼出了什麼問題以及它發生了什麼。

如果有人能指出錯誤或任何其他建議，那將是非常有用的。

這裏是寫在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應用於列時，上述操作纔有效。如果輸入是列，size(signal,1)肯定會給你輸入矢量的長度。但是，如果您的輸入是行，那麼size(signal,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 

---

正如你所看到的，結果相匹配。在我看來很好！

---

只是爲了確保我們是一致的，這是完整的代碼爲您的功能列表，以及我所做的修改：

% 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