离散余弦变换DCT

大家好，欢迎来到IT知识分享网。

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

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

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

基本介绍

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

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

特点总结：

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

数学公式

一维DCT：

二维DCT：

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 

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

二维DCT反变换：

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反变换恢复图');码片 

例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) 

在 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) 

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

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