宝塔服务器面板,一键全能部署及管理,送你10850元礼包,点我领取

背景与原理
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)=NM1x=0N1y=0M1f(x,y)eN2πiuxeM2πivyf(x,y)=NM1u=0N1v=0M1f(u,v)eN2πiuxeM2πivy
傅里叶变换包含复数运算,其运算复杂度和存储长度都超过实数运算。为了简化上述过程,同时达到更好的变换效果,余弦变换应运而生。

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

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

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

假如序列中不包含零点,自然按照以下方式延拓:
f(1,0)=f(−1,0),f(0,1)=f(0,−1),对称中心为(0,0)f(1,0)=f(-1,0),f(0,1)=f(0,-1),对称中心为(0,0)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)f(0,0)=f(-1,0),f(0,0)=f(0,-1),对称中心为(-1/2,-1/2)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)=N1x=0N1y=0N1f(x,y)cos[Nπu(x+21)]cos[Nπv(y+21)]f(u,v)=N1x=0N1y=0N1F(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)=2N1x=0N1y=0N1f(x,y)cos[Nπu(x+21)]cos[Nπv(y+21)]f(u,v)=2N1x=0N1y=0N1F(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}eN2πiuxeM2πivye2πiNuxe2πiMvye^{\frac{2\pi i}{N}ux}e^{\frac{2\pi i}{M}vy}eN2πiuxeM2πivy相差一个负号;相应的,在余弦变换中,正逆变换的变换核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)cos(x)=cos(-x)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)=N1x=0N1y=0N1f(x,y)cos[2Nπu(2x+1)]cos[2Nπv(2y+1)]
也就是将12\frac{1}{2}21 乘入了整个分式中。虽然这样写只是变了一种形式,不影响结果,但破坏了原本的几何意义(或者说偶延拓过程),谁知道这个 2x+12x+12x+1 对应着什么东西。

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

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

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

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

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

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

二、DCT在硬件层面的实现
我们先不谈DSP等专用芯片,单说CPUGPU

计算机的运算可以分为两类,整数运算浮点运算;浮点运算又可以分为两类,单精度浮点运算双精度浮点运算

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)=161x=015y=015cos[8πu(x+21)]cos[8πv(y+21)]f(x,y)=161u=015v=015cos[8πu(x+21)]cos[8πv(y+21)]在整个编码过程中,离散余弦变换的作用如下:

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

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


以上博文转载自:

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