离散余弦变换DCT

离散余弦变换（Discrete Cosine Transform, DCT）

百度百科的解释：
离散余弦变换(DCT)是与傅里叶变换相关的一种变换，它类似于离散傅里叶变换(DFT)。但是只使用实数，DCT相当于一个长度大概是DFT两倍的离DFT，这个离散傅里叶变换是对一个实偶函数进行的（因为一个实偶函数的傅里叶变换仍然是一个实偶函数），在有些变形里面需要将输入或者输出的位置移动半个单位(DCT有8种标准类型，其中4种是常见的)。

离散傅里叶变换需要进行复数运算，尽管有FFT可以提高运算速度，但在图像编码、特别是在实时处理中非常不便。离散傅里叶变换在实际的图像通信系统中很少使用，但它具有理论的指导意义。根据离散傅里叶变换的性质，实偶函数的傅里叶变换只含实的余弦项，因此构造了一种实数域的变换——离散余弦变换(DCT)。

基本介绍

通过研究发现，DCT除了具有一般的正交变换性质外，其变换阵的基向量很近似于Toeplitz矩阵的特征向量，后者体现了人类的语言、图像信号的相关特性。因此，在对语音、图像信号变换的确定的变换矩阵正交变换中，DCT变换被认为是一种准最佳变换。在近年颁布的一系列视频压缩编码的国际标准建议中，都把 DCT 作为其中的一个基本处理模块。

DCT变换较DFT变换具有更好的频域能量聚集度，对于不重要的频域区域和系数就能够直接裁剪掉，因此，DCT变换非常适合于图像压缩算法的处理，例如现在大名鼎鼎的jpeg就是使用了DCT作为图像压缩算法。

DCT除了上述介绍的几条特点，即：实数变换、确定的变换矩阵、准最佳变换性能外，二维DCT还是一种可分离的变换，可以用两次一维变换得到二维变换结果。

特点总结：

• DCT是一种有损变换，但特点是其避免了傅立叶变换中的复数实部虚部问题。
• DCT使信号低频区集中在左上角便于裁剪获取低频信息，而自然信号中的能量多集中在低频区(即属性变化平缓的区域，如一张图像中像素值变化平缓的区域往往包含主要的色彩信息而变化剧烈的区域往往是线条边界和色彩交替区)，故常用于信号和图像处理中数据的有损压缩。

数学公式

一维DCT：

一维DCT变换共有8种形式，其中最常用的是第二种形式，由于其运算简单、适用范围广。我们在这里只讨论这种形式，其表达式如下：

其中，f(i)为原始的信号，F(u)是DCT变换后的系数，N为原始信号的点数，c(u)可以认为是一个补偿系数，可以使DCT变换矩阵为正交矩阵。

二维DCT：

二维DCT变换其实是在一维DCT变换的基础上再做了一次DCT变换，其公式如下：

由公式我们可以看出，上面只讨论了二维图像数据为方阵的情况，在实际应用中，如果不是方阵的数据一般都是补齐之后再做变换的，重构之后可以去掉补齐的部分，得到原始的图像信息。
另外，由于DCT变换高度的对称性，在使用Matlab进行相关的运算时，我们可以使用更简单的矩阵处理方式：

DCT 系数F(u,v)可视为应用于每个基函数的权重。对于 8×8 矩阵，该图可表示为 64 个基函数。

水平频率从左到右递增，垂直频率从上到下递增。左上角的常量值基函数通常称为 DC 基函数，对应的 DCT 系数 B00 通常称为 DC 系数。

DCT变换与IDCT变换，MATLAB代码实现：

clear;
clc;

% 正变换
X=round(rand(4)*100)   %产生随机矩阵
A=zeros(4);
for i=0:3
    for j=0:3
        if i==0
            a=sqrt(1/4);
         else
             a=sqrt(2/4);
         end            
         A(i+1,j+1)=a*cos(pi*(j+0.5)*i/4);
    end
end

Y=A*X*A'        %DCT变换

%反变换
for i=0:3
    for j=0:3
         if i==0
             a=sqrt(1/4);
         else
             a=sqrt(2/4);
         end            
         A(i+1,j+1)=a*cos(pi*(j+0.5)*i/4); %生成变换矩阵
     end
end

X1=A'*Y*A  %DCT反变换恢复的矩阵

% Matlab版
YY=dct2(X)      %Matlab自带的dct变换
XX=idct2(YY)    %Matlab自带的idct逆变换

运行结果为：

 1 X =
 2 
 3     42    66    68    66
 4     92     4    76    17
 5     79    85    74    71
 6     96    93    39     3
 7 
 8 
 9 Y =
10 
11   242.7500   48.4317   -9.7500   23.5052
12   -12.6428  -54.0659    7.4278   22.7950
13    -6.2500   10.7158  -19.7500  -38.8046
14    40.6852  -38.7050  -11.4653  -45.9341
15 
16 
17 YY =
18 
19   242.7500   48.4317   -9.7500   23.5052
20   -12.6428  -54.0659    7.4278   22.7950
21    -6.2500   10.7158  -19.7500  -38.8046
22    40.6852  -38.7050  -11.4653  -45.9341

由上面的结果我们可以看出，我们采用的公式的方法和Matlab自带的dct变化方法结果是一致的，所以验证了我们方法的正确性。
如果原始信号是图像等相关性较大的数据的时候，我们可以发现在变换之后，系数较大的集中在左上角，而右下角的几乎都是0，其中左上角的是低频分量，右下角的是高频分量，低频系数体现的是图像中目标的轮廓和灰度分布特性，高频系数体现的是目标形状的细节信息。DCT变换之后，能量主要集中在低频分量处，这也是DCT变换去相关性的一个体现。

之后在量化和编码阶段，我们可以采用“Z”字形编码，这样就可以得到大量的连续的0，这大大简化了编码的过程。

二维DCT反变换：

在图像的接收端，根据DCT变化的可逆性，我们可以通过DCT反变换恢复出原始的图像信息，其公式如下：

同样的道理，我们利用之前的矩阵运算公司可以推导出DCT反变换相应的矩阵形式：

Matlab代码：

clear;
clc;
X=[  61    19    50    20
     82    26    61    45
     89    90    82    43
     93    59    53    97] %原始的数据

A=zeros(4);
for i=0:3
     for j=0:3
         if i==0
             a=sqrt(1/4);
         else
             a=sqrt(2/4);
         end            
         A(i+1,j+1)=a*cos(pi*(j+0.5)*i/4); %生成变换矩阵
     end
 end
 Y=A*X*A'   %DCT变换后的矩阵
 X1=A'*Y*A  %DCT反变换恢复的矩阵

运行结果：

 1 X =
 2 
 3     61    19    50    20
 4     82    26    61    45
 5     89    90    82    43
 6     93    59    53    97
 7 
 8 
 9 Y =
10 
11   242.5000   32.1613   22.5000   33.2212
12   -61.8263    7.9246  -10.7344   30.6881
13   -16.5000  -14.7549   22.5000   -6.8770
14     8.8322   16.6881  -35.0610   -6.9246
15 
16 
17 X1 =
18 
19    61.0000   19.0000   50.0000   20.0000
20    82.0000   26.0000   61.0000   45.0000
21    89.0000   90.0000   82.0000   43.0000
22    93.0000   59.0000   53.0000   97.0000

我们可以看到反变换后无损的恢复了原始信息，所以证明了方法的正确性。但是在实际过程中，需要量化编码或者直接舍弃高频分量等处理，所以会出现一定程度的误差，这个是不可避免的。

应用举例

使用离散余弦变换的图像压缩
例1：
在实际的图像处理中，DCT变换的复杂度其实是比较高的，所以通常的做法是，将图像进行分块，然后在每一块中对图像进行DCT变换和反变换，在合并分块，从而提升变换的效率。具体的分块过程中，随着子块的变大，算法复杂度急速上升，但是采用较大的分块会明显减少图像分块效应，所以，这里面需要做一个折中，在通常使用时，大都采用的是8*8的分块。
Matlab的 blkproc 函数可以帮我们很方便的进行分块处理，下面给出我们的处理过程：

 1 clear;
 2 clc;
 3  
 4 X=imread('pepper.bmp');
 5 X=double(X);
 6 [a,b]=size(X);
 7 Y=blkproc(X,[8 8],'dct2');
 8 X1=blkproc(Y,[8 8],'idct2');
 9  
10 figure
11 subplot(1,3,1);
12 imshow(uint8(X));
13 title('原始图');
14  
15 subplot(1,3,2);
16 imshow(uint8(Y));
17 title('分块DCT变换图');
18  
19 subplot(1,3,3);
20 imshow(uint8(X1));
21 title('分块DCT恢复图');
22  
23 Y1=dct2(X);
24 X10=idct2(Y1);
25  
26 figure
27 subplot(1,3,1);
28 imshow(uint8(X));
29 title('原始图');
30  
31 subplot(1,3,2);
32 imshow(uint8(Y1));
33 title('DCT变换图');
34  
35 subplot(1,3,3);
36 imshow(uint8(X10));
37 title('DCT反变换恢复图');码片

结果：


从图中，我们可以明显看出DCT变换与分块DCT变换在使用时的区别。

例2：

因为噪声主要存在于高频信息中，对高频信息进行适当抑制，可以起到图像去噪的作用，这里采用简单高频抑制方法，可以降噪但也会丢失细节，中间处理的方法还有很多就不一一列举，MATLAB代码如下：

%读取图像
X=imread('lena.jpg'); 
X=rgb2gray(X);
%读取图像尺寸
[m,n]=size(X); 
%给图像加噪
Xnoised=imnoise(X,'gaussian',0.01); 
%输出加噪图像
subplot(121); 
imshow(Xnoised);
%DCT变换
Y=dct2(Xnoised); 
I=zeros(m,n);
%高频抑制
I(1:m/3,1:n/3)=1; 
Ydct=Y.*I;
%逆DCT变换
Y=uint8(idct2(Ydct)); 
%结果输出
subplot(122);
imshow(Y);

Python代码如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# DCT hyoer-parameter 超参数
T = 8
K = 8
channel = 3

# DCT weight
def w(x, y, u, v):
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))

# DCT
def dct(img):
    H, W, _ = img.shape

    F = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for v in range(T):
                    for u in range(T):
                        for y in range(T):
                            for x in range(T):
                                F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)

    return F


# IDCT
def idct(F):
    H, W, _ = F.shape

    out = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for y in range(T):
                    for x in range(T):
                        for v in range(K):
                            for u in range(K):
                                out[y+yi, x+xi, c] += F[v+yi, u+xi, c] * w(x,y,u,v)

    out = np.clip(out, 0, 255)
    out = np.round(out).astype(np.uint8)

    return out



# Read image
img = cv2.imread("imori.jpg").astype(np.float32)

# DCT
F = dct(img)

# IDCT
out = idct(F)

# Save result
cv2.imshow("result", out)
cv2.waitKey(0)
cv2.imwrite("out.jpg", out)

例3：DCT
此示例说明如何使用离散余弦变换 (DCT) 压缩图像。该示例计算输入图像中 8×8 个数据块的二维 DCT，丢弃（设置为零）每个数据块中 64 个 DCT 系数中除 10 个以外的所有系数，然后使用每个数据块的二维逆 DCT 重构图像。该示例使用变换矩阵计算方法。

在 JPEG 图像压缩算法中使用 DCT。将输入图像分成 8×8 或 16×16 个数据块，并对每个数据块计算二维 DCT。然后对 DCT 系数进行量化、编码和传输。JPEG 接收机（或 JPEG 文件读取器）对量化的 DCT 系数进行解码，计算每个数据块的逆二维 DCT，然后将这些数据块一起放回单个图像中。对于典型图像，许多 DCT 系数的值接近于零。可以丢弃这些系数，而不会严重影响重构图像的质量。

I = imread('cameraman.tif'); %将图像读入工作区
I = im2double(I); %转换为 double 类

%计算图像中 8×8 个数据块的二维 DCT。函数 dctmtx 返回 N×N DCT 变换矩阵。
T = dctmtx(8);
dct = @(block_struct) T * block_struct.data * T';
B = blockproc(I,[8 8],dct);

%丢弃每个数据块中 64 个 DCT 系数的大部分系数，仅保留 10 个。
mask = [1   1   1   1   0   0   0   0
        1   1   1   0   0   0   0   0
        1   1   0   0   0   0   0   0
        1   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0];
B2 = blockproc(B,[8 8],@(block_struct) mask .* block_struct.data);

%使用每个数据块的二维逆 DCT 重构图像。
invdct = @(block_struct) T' * block_struct.data * T;
I2 = blockproc(B2,[8 8],invdct);

imshow(I)
figure
imshow(I2)

结果：并排显示原始图像和重构图像。尽管几乎 85% 的 DCT 系数被丢弃，导致重构图像的质量有所下降，但它仍是清晰可辨的。

注：
使用 Image Processing Toolbox™ 软件有两种计算 DCT 的方法。

• dct2 函数：使用基于 FFT 的算法对大型输入实现快速计算。
• dctmtx函数： DCT 变换矩阵，该矩阵由函数 dctmtx 返回，对于小型方阵输入（如 8×8 或 16×16）可能更高效。

M×M 变换矩阵 T 由下式给出：

对于 M×M 矩阵 A，TA 是 M×M 矩阵，其列包含由 A 的列组成的一维 DCT。A 的二维 DCT 可以计算为 B=TA*T’。由于 T 是实数正交矩阵，它的逆矩阵与它的转置矩阵相同。因此，B 的逆二维 DCT 由 T’BT 给出。