离散余弦变换原理及实现过程【转载】

背景与原理
1974年，K. R. Rao、N. Ahmed、T. Natarajan三位教授创立了离散余弦变换（Discrete Cosine Transform, DCT）。在数字信号、数字图像处理领域，离散余弦变换的效果能够接近理论上的最佳变换——Kahunen-Loeve变换（K-L变换）。

以下将介绍DCT的相关背景，并从算法、硬件、应用三个层面进行概述。

1807年，法国数学家、物理学家傅里叶（Jean Baptiste Joseph Fourier）提出了傅里叶变换（Fourier Transform, FT）。傅里叶变换的形式有很多种，归一化的二维离散傅里叶变换（Discrete Fourier transform, DFT）可以写成如下形式：
F(u,v)=1NM∑x=0N−1∑y=0M−1f(x,y)e−2πiNuxe−2πiMvyf(x,y)=1NM∑u=0N−1∑v=0M−1f(u,v)e2πiNuxe2πiMvyF(u,v)=\frac{1}{\sqrt{NM}}\sum^{N-1}_{x=0}\sum^{M-1}_{y=0}f(x,y)e^{-\frac{2\pi i}{N}ux}e^{-\frac{2\pi i}{M}vy}\\ f(x,y)=\frac{1}{\sqrt{NM}}\sum^{N-1}_{u=0}\sum^{M-1}_{v=0}f(u,v)e^{\frac{2\pi i}{N}ux}e^{\frac{2\pi i}{M}vy} F(u,v)=NM​1​x=0∑N−1​y=0∑M−1​f(x,y)e−N2πi​uxe−M2πi​vyf(x,y)=NM​1​u=0∑N−1​v=0∑M−1​f(u,v)eN2πi​uxeM2πi​vy
傅里叶变换包含复数运算，其运算复杂度和存储长度都超过实数运算。为了简化上述过程，同时达到更好的变换效果，余弦变换应运而生。

从傅里叶变换到离散余弦变换，需要一些数学理论的支持。在给定区间内满足狄利赫里条件的连续实对称函数，可以展开成仅含有余弦项的傅里叶级数。

对于定义在正实数域上的函数，可以通过偶延拓或奇延拓，满足上述条件。但如果函数的定义域包含零点，情况则稍有些复杂。

以一个二维离散函数$f(x,y)(x,y=0,1,...,N-1)$为例，对其进行偶延拓。

假如序列中不包含零点，自然按照以下方式延拓：
$f(1,0)=f(-1,0),f(0,1)=f(0,-1),对称中心为(0,0)$

由于序列中包括零点，考虑零点后的延拓方式：
$f(0,0)=f(-1,0),f(0,0)=f(0,-1),对称中心为(-1/2,-1/2)$

因此，按照上述方法延拓后，归一化的二维离散余弦变换可以写成如下形式：
when(x,y)or(u,v)=(0,0)F(u,v)=1N∑x=0N−1∑y=0N−1f(x,y)cos[πNu(x+12)]cos[πNv(y+12)]f(u,v)=1N∑x=0N−1∑y=0N−1F(u,v)cos[πNu(x+12)]cos[πNv(y+12)]when (x,y) or (u,v)=(0,0)\\ F(u,v)=\frac{1}{N}\sum^{N-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)cos[\frac{\pi}{N}u(x+\frac{1}{2})]cos[\frac{\pi}{N}v(y+\frac{1}{2})]\\ f(u,v)=\frac{1}{N}\sum^{N-1}_{x=0}\sum^{N-1}_{y=0}F(u,v)cos[\frac{\pi}{N}u(x+\frac{1}{2})]cos[\frac{\pi}{N}v(y+\frac{1}{2})] when(x,y)or(u,v)=(0,0)F(u,v)=N1​x=0∑N−1​y=0∑N−1​f(x,y)cos[Nπ​u(x+21​)]cos[Nπ​v(y+21​)]f(u,v)=N1​x=0∑N−1​y=0∑N−1​F(u,v)cos[Nπ​u(x+21​)]cos[Nπ​v(y+21​)]

when(x,y)or(u,v)≠(0,0)F(u,v)=12N∑x=0N−1∑y=0N−1f(x,y)cos[πNu(x+12)]cos[πNv(y+12)]f(u,v)=12N∑x=0N−1∑y=0N−1F(u,v)cos[πNu(x+12)]cos[πNv(y+12)]when (x,y) or (u,v)\neq(0,0)\\ F(u,v)=\frac{1}{2N}\sum^{N-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)cos[\frac{\pi}{N}u(x+\frac{1}{2})]cos[\frac{\pi}{N}v(y+\frac{1}{2})]\\ f(u,v)=\frac{1}{2N}\sum^{N-1}_{x=0}\sum^{N-1}_{y=0}F(u,v)cos[\frac{\pi}{N}u(x+\frac{1}{2})]cos[\frac{\pi}{N}v(y+\frac{1}{2})] when(x,y)or(u,v)​=(0,0)F(u,v)=2N1​x=0∑N−1​y=0∑N−1​f(x,y)cos[Nπ​u(x+21​)]cos[Nπ​v(y+21​)]f(u,v)=2N1​x=0∑N−1​y=0∑N−1​F(u,v)cos[Nπ​u(x+21​)]cos[Nπ​v(y+21​)]

在傅里叶变换中，正逆变换的变换核e−2πiNuxe−2πiMvye^{-\frac{2\pi i}{N}ux}e^{-\frac{2\pi i}{M}vy}e−N2πi​uxe−M2πi​vy和e2πiNuxe2πiMvye^{\frac{2\pi i}{N}ux}e^{\frac{2\pi i}{M}vy}eN2πi​uxeM2πi​vy相差一个负号；相应的，在余弦变换中，正逆变换的变换核cos[πNu(x+12)]cos[πNv(y+12)]cos[\frac{\pi}{N}u(x+\frac{1}{2})]cos[\frac{\pi}{N}v(y+\frac{1}{2})]cos[Nπ​u(x+21​)]cos[Nπ​v(y+21​)]也应相差一个负号。由于 $cos(x)=cos(-x)$ ，所以余弦变换的正逆变换在形式上具有一致性。

顺带一提，在某些教材上，将上述形式的离散余弦变换写成类似于下面的样子：
F(u,v)=1N∑x=0N−1∑y=0N−1f(x,y)cos[πu(2x+1)2N]cos[πv(2y+1)2N]F(u,v)=\frac{1}{N}\sum^{N-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)cos[\frac{\pi u(2x+1)}{2N}]cos[\frac{\pi v(2y+1)}{2N}] F(u,v)=N1​x=0∑N−1​y=0∑N−1​f(x,y)cos[2Nπu(2x+1)​]cos[2Nπv(2y+1)​]
也就是将12\frac{1}{2}21​ 乘入了整个分式中。虽然这样写只是变了一种形式，不影响结果，但破坏了原本的几何意义（或者说偶延拓过程），谁知道这个 $2x+1$ 对应着什么东西。

一、DCT在算法层面的实现
为了简化离散余弦变换的计算过程，需要算法上进行优化。优化可以从两方面入手，一是降低计算次数，一是将浮点运算转化为整数运算。

在计算次数方面，对于$N$次离散余弦变换，需要o(N2)o(N^2)o(N2)次运算，通过快速算法，能够将运算降低至o(Nlog2N)o(Nlog_2N)o(Nlog2​N) 次。

在整数运算方面，通过量化实现整数近似，并且在量化的同时保证变换的可逆性（不然变过去之后，就变不回来了），将离散余弦变换过程中的浮点运算转化为整数运算，大致过程如下：

为了熟悉和理解上的方便，将离散余弦变换表达式写成矩阵形式：
Y=CXCTY=CXC^T Y=CXCT所谓整数近似，并不是单纯地将浮点数四舍五入或者直接截断得到整数（这样做后，再进行逆变换就变不回来了），而是将浮点运算归入量化过程（在逆变换时进行反量化，保证与原始数据的一致性），让离散余弦变换过程仅包含整数运算。

这样做看似是朝三暮四，但在实际应用中，量化和量化中的浮点运算是必不可少的。反正总归是要量化，这样可以将两次浮点预算合并为一次，对优化算法仍旧是有意义的。

进行量化的整合后，离散余弦变换表达式为：
Y=(CfXCfT)⨀EfY=(C_fXC^T_f)\bigodot E_f Y=(Cf​XCfT​)⨀Ef​其中，$⨀$表示矩阵的点乘。完成上述整数化后，不仅将浮点运算归入量化过程，而且可以并行计算，加快速度。

二、DCT在硬件层面的实现
我们先不谈DSP等专用芯片，单说CPU与GPU。

计算机的运算可以分为两类，整数运算和浮点运算；浮点运算又可以分为两类，单精度浮点运算和双精度浮点运算。

GPU拥有大量并行流处理器，但缺少串行逻辑控制机构，擅长进行浮点运算，对于更注重溢出检查的整数运算，反而是GPU的弱项；CPU作为中央处理器，则是两种运算通吃。

进一步考虑浮点运算，单精度浮点运算的复杂度低于双精度浮点运算。市面上常见的游戏显卡和TaiTan X，均是注重单精度浮点运算，而限制了双精度浮点运算。开普勒架构的游戏显卡，双精度浮点运算能力是单精度浮点运算的1/24；同时代的CPU，这一比值大约是1/2。

一般而言，对于物理建模与模拟（比如流体力学模拟、量子化学计算）、3D建模（其实也是一种意义上的物理建模），需要双精度浮点运算；对于一般的图形渲染（包括游戏渲染、视频渲染）、机器学习，需要单精度浮点运算。所以游戏显卡削弱双精度浮点运算，TaiTan X一个劲提升单精度浮点运算，以及1080Ti被视为平民机器学习神器，都是有其道理的。

回到离散余弦变换。的确，GPU的并行运算速度超过CPU。经过优化的整数离散余弦变换算法利于并行计算，更是突出了GPU的并行优势，但GPU并不擅长整数运算。此外，在实际应用中，离散余弦变换过程伴随着串行逻辑，这也GPU不能胜任的。随着GPU计算的火热，如何让GPU完成离散余弦变换成了学界热点。但目前来看，让CPU扮演主要角色，把一部分运算交给GPU实现加速，应当是比较现实的考虑。

三、DCT在应用层面的实现
以视频编码为例。1984年，H.261编码标准开始研究。1990年，H.261编码标准由国际电信联盟（ITU-T）发布。1993年，MPEG-1编码标准由国际标准化组织/国际电子学委员会（ISO/IEC）发布。自此，H.26X系列与MPEG-X系列开始了各自的更新换代，当然期间也有合作交易，比如大家可能比较熟悉的MPEG-2与大家可能不熟悉的H.262，其实是一个东西，又比如MPEG-4，包括了（早期的）MPEG-4、MPEG-4/AVC（=H.264）、MPEG-4/HEVC（=H.265）三个标准。

诶，好像就差H.263与MPEG-3没有提。前者随着技术的进步退出了历史舞台，后者则被毙掉了…

从H.261开始，离散余弦变换就被应用于视频编码中。当技术发展到H.264与H.265，就具体的数学计算而言，H.261创立的分块编码被后续标准沿用，H.264在亮度平面上的基本分块为8×8像素块，于是离散余弦变换是这个样子（代入N=8）：
F(u,v)=116∑x=015∑y=015cos[π8u(x+12)]cos[π8v(y+12)]f(x,y)=116∑u=015∑v=015cos[π8u(x+12)]cos[π8v(y+12)]F(u,v)=\frac{1}{16}\sum^{15}_{x=0}\sum^{15}_{y=0}cos[\frac{\pi}{8}u(x+\frac{1}{2})]cos[\frac{\pi}{8}v(y+\frac{1}{2})]\\ f(x,y)=\frac{1}{16}\sum^{15}_{u=0}\sum^{15}_{v=0}cos[\frac{\pi}{8}u(x+\frac{1}{2})]cos[\frac{\pi}{8}v(y+\frac{1}{2})] F(u,v)=161​x=0∑15​y=0∑15​cos[8π​u(x+21​)]cos[8π​v(y+21​)]f(x,y)=161​u=0∑15​v=0∑15​cos[8π​u(x+21​)]cos[8π​v(y+21​)]在整个编码过程中，离散余弦变换的作用如下：

首先，按照图像处理的一般流程，进行采样、整量，完成连续数据到离散数据的量化。根据前面提到的算法优化，量化过程整合了离散余弦变换的浮点运算。

接下来，进行整数离散余弦变换。做这一变换有什么用？一方面，从图像处理的整体流程而言，变换后便于后续处理；另一方面，从编码的角度而言，变换后使图像信息集中，在数学上体现为描述关键信息的系数变少，相应的，所需存储空间降低，达到降低视频体积的目的。

---

以上博文转载自：

1. https://zhuanlan.zhihu.com/p/33845296