基于CUDA技术的DCT并行算法研究与实现

张金霜 黄旭彬

摘要：JPEG有损压缩算法以DCT变换算法为核心，但DCT变换算法较为耗时，为提升图像压缩效率，提出利用基于GPU平台的CUDA技术对DCT算法做并行优化。通过分析DCT变换算法的原理，发现DCT算法具有很好的并行性，优化空间较大，于是利用CUDA技术实现高速DCT并行算法。实验结果表明，在一张2592x2592分辨率的图像做DCT变换，加速比能达到124.5，算法效率得到明显提升，且压缩效果无显著差异。

关键词：CUDA;GPU;并行程序设计;DCT算法;图像压缩

中图分类号：TP393 文献标识码：A

文章编号：1009-3044（2020）33-0198-04

开放科学（资源服务）标识码（OSID）：

1 背景

离散余弦变换（ Discrete Cosine Transform，DCT）类似于傅立叶变换，由Ahmed在1974年首次提出。离散余弦变换经常被用作信号和图像处理，如对信号和图像进行有损数据压缩。在JPEG图像压缩算法中，有一种是以离散余弦变换为基础的有损压缩算法[1]。作为JPEG压缩算法的核心，DCT是其中最耗时的部分，且占有较大的工作量比重。

目前GPU作为一种擅长处理大规模并行数据计算的设备，已受到越来越多人的关注。相较于CPU，GPU有着强大的计算能力和超高度存储器带宽，可以使用数量较多的执行单元，实现数据间的并发计算。而CUDA技术能充分利用GPU的并行计算能力，且降低了开发人员的使用门槛，基于GPU的CUDA并行加速技术也广泛应用于图像处理领域[2-4]。因此，本文将DCT串行算法改为基于CPU/GPU的异构算法，通过GPU内部大量独立的计算单元，实现计算加速。

2 CUDA体系结构

CUDA是由显卡厂商NVIDIA推出的通用并行计算框架，该架构能使有效调用GPU资源解决复杂的计算问题。CUDA允许开发人员使用C/C++语言和CUDA扩展库的形式编写程序，避免了传统方法中将通用计算转换为图形处理的过程，大大降低程序设计的难度和系统维护成本。

CUDA实现了CPU和GPU的混合架构编程。在CUDA架构下，程序被分成两个部分：host端程序和device端程序。host端程序运行在CPU上;device端程序又叫“kernel”，运行在GPU上。通常先由host端程序准备数据，对设备进行必要的初始化，把数据从内存传送到显存中;执行device端程序，待GPU运行完后，host端程序再从显存中把结果取回内存中，过程如图1所示[5]。

CUDA的线程模型分三个层次：thread、block、grid。其中，kernel以线程网格（grid）的形式组织，每个线程网格由若干个线程块（block）组成，而每个线程块又由若干个线程（thread）组成。

block内的thread之间是并行执行的，grid内的block之間也是并行执行的，而grid与grid之间是串行执行的，即GPU同时只能处理一个grid。grid包含一组block，grid是kernel在GPU上的一次执行，也就是说，一个kernel对应着一个grid。grid内部的block可以按照一维或者二维来组织，block之间共享全局内存。

3 DCT算法理论

二维频率与矩阵c（u，v）的左上角代表二维正交展开中的低阶项，它代表图像的低频分量;c（u，v）的右下角代表二维正交展开中的高阶项，它代表图像的高频分量。图像的低频分量，反映图像灰度的缓慢变化;相反，图像的高频分量反映图像灰度的跳变，又称轮廓或边缘[6]。

经过DCT变换后的矩阵数据自然数为频率系数，这些系数中，位于左上角的F（0，0）的值最大，称为直流系数（DC coeffi-cient）;其余的63个频率系数多为一些接近于0的正负浮点数，称为交流系数（AC coeffcients）[7]。利用这一特点就可以实现图像数据压缩，对其进行量化，仅仅保留低频分量的左上角，其余分量抛弃，这样就可以达到图像数据压缩的目的。在反变换的时候，抛弃的点用零填充。

4 基于CUDA实现DCT算法加速

4.1可行性分析

CUDA架构的优势在于处理数据并行的问题。通过生成尽可能多的线程来处理，每一个线程的处理是基本相同的，而处理的数据是可以并行的。符合这种特点的算法，能在CUDA架构下最大化并行执行。

根据JPEG标准，图像被分成一个个8x8的像素块进行DCT变换。针对每一个8x8的像素块，它们的处理是相同的，都是进行DCT变换。而每一个8x8的像素块是相互独立的，块与块之间没有关联。而变换矩阵是确定的，与变换对象无关。因此对各8x8像素块进行DCT变换是数据并行的。

4.2 算法总体流程

1）算法从host端开始，首先在内存中分配两个具有与被处理图像像素数据相同大小的空间，分别用于存放被处理图像的源数据和处理结果，并把图像像素数据读人到内存中。图像要在GPU中进行处理，于是在显存中申请两个具有与被处理图像像素数据相同大小的空间，分别用于存放图像源数据和处理结果，并一次性把图像像素数据从内存复制到显存中。

2）把每一个8x8块读人到共享存储器中。采用多线程读人，每个线程读人一个像素。读人完成后，进行DCT变换的计算。采用多线程进行计算，每个线程计算一个系数。

3）DCT变换完成后，把结果写入到全局存储器中，同样是每个线程写入一个系数。至此device端的工作完成。

4）host端再把结果数据从显存复制到内存中，DCT变换算法完成。

4.3 算法详细设计

1）线程划分

把整个图像作为一个grid，由于对图像采取的是8x8块的DCT变换，块与块之间可以独立进行处理，不同块没有关联。因此，这里采取每个8x8块用一个block来处理，而block中包含8x8个thread（共64个thread），每个thread处理一个像素的计算。工作过程如图2所示。

2）存储器分配

根据各存储器的特点，先把图像数据载人到全局存储器中。由于全局存储器的访问速度慢，而且没有缓冲，每一个block内，各thread要取的数据有重复之处，如果线程直接从全局存储器取数据，将会遭遇性能瓶颈。因此，由各thread从全局存储器中读人到共享存储器中，各自线程读取自身坐标对应的像素值。在各thread读入像素值后，要进行栅障同步，以保证同一个block内所有像素都已经读入，这样才能进行下一步的DCT变换。

为了进一步优化，还可以把全局存储器的中数据映射为纹理，纹理有片内缓冲，而且能进行二维寻址，能简化DCT算法的寻址，提高访存的速度。同时，纹理数据按16的倍数对齐，以符合“合并访问”的要求，进一步提高访存效率。

此外，DCT变换所需的变换矩阵A要先载人到常数存储器，常数存储器同样具有片内缓冲。在代码中，用一constant一声明限定符直接声明为常数存储器变量即可[8]。

3）并行DCT算法实现

根据等式（5），DCT变换要计算两次矩阵乘法，分别是ATXArX和（ArX）A。矩阵相乘的算法，可以通过线程的索引来确定其计算的数据地址。在CUDA中，线程是能获得自身在block中的索引的，基于这一特点，在矩阵乘法运算中，可以很方便地在共享存储器中寻址。每个矩阵都需要用到转置矩阵A，而矩阵A被载入到常数存储器，GPU能对常数存储器进行片内缓冲，访问效率很高。

以下为DCT变换中两次求矩阵乘积的核心代码：

//计算ATX，结果保存到CurBlockLoca12

float curelem=0;//临时保存当前元素的值

int DCTv8matrixlndex=ty;//转置矩阵A1索引

int CurBlockLocalllndex=tx;//x矩阵的索引

for （int i=0; i

{

curelem+= DCTv8matrix[DCTv8matrixlndex] * Cur-BlockLocaII[CurBlockLocalllndexl;

DCTv8matrixlndex+=BLOCK_SIZE; //A矩阵的列为AT矩阵的行

CurBlockLocalllndex+=BLO CK_SIZE;

）

CurBlockLoca12[ （ty << BLOCK—SIZE—.LOG2） + tx]=curelem;/，保存结果

_syncthreads0;//等待矩阵的所有元素计算完毕

//计算（ATX） A（ATX）A，结果保存到CurBlockLocall

curelem=0;//临时保存当前元素的值

int CurBlockLoca121ndex

=

（ty

*BLOCK_SIZE）;//（ATX ）（ATX）矩阵索引

DCTv8matrixlndex=tx;//A矩阵索引

for （int i=0; i

（

curelem += CurBlockLoca12[CurBlockLoca121ndex] *DCTv8matrix[DCTv 8matrixlndexl;

CurBlockLocal21ndex+=1;//取行，因此步长为1

DCTv8matrixlndex+=BLOCK_SIZE;，/取列，因此步长为BLOCK__ SIZE

}

CurBlockLocall[（ty<

_syncthreads0;//栅障同步，确保所有系数计算完成

5 实验结果分析

5.1 数据分析

为了测试应用CUDA技术后的运行效率，本文编写了两个版本的程序，分别在CPU运行（单线程），和在GPU上运行（基于CUDA技术）。通过分别调用两个程序，针对同一幅图片进行DCT变换，并计算其执行时间。

在计算执行时间上，本文采用CUDA的cutStartTimer（tim-er）、cutStopTimer（timer）和cutGetTimerValue（timer）函数，分别表示开始计时、停止计时和读取时间，计时精度小于Ims。

硬件上，CPU采用Intel Core i5-7200，主频为2.5GHz;GPU采用NVIDIA公司的GeForce GTX660，支持CUDA 4.0以上驱动版本，显存容量2G。

经过测试，DCT算法在CPU和GPU上执行的效率对比如表1所示。从表中数据可以看出，应用CUDA技术加速DCT变换，加速比最高能达到124.5。

5.2 程序演示

本文通過模块化的方式，把DCT变换算法这个功能作为一个独立模块。DCT变换算法的CUDA实现版本和CPU实现版本（C语言洛编译成一个DLL文件。主程序模块则由C#语言编写，编译成一个可执行文件。

以下程序分别演示了应用CUDA技术在GPU上进行图像压缩，和同样的算法在CPU上进行计算的效果。左边的图是原图，中间的图是经过DCT变换后的系数，右边的图是经过DCT——量化——逆量化-IDCT后的结果。

6 结束语

本文重点分析了DCT变换算法的原理，得出DCT变换算法的矩阵式表示，并成功利用CUDA技术实现了高速DCT变换算法。实验表明，相较于CPU上的串行DCT变换算法，本文算法的执行效率得到了很大提升。

本文也提供了一种程序优化思路——对于耗时的串行程序，只要满足并行性的条件，都可尝试采用CUDA技术进行改进。而把CUDA部分的程序编译成DLL文件，可灵活地被第三方项目调用，使得CUDA技术适用于更多技术场景。

参考文献：

[1]万丰，苑豪杰，宫威，基于离散余弦变换的数字图像压缩技术研究[J].自动化应用，2020（3）：65-67.

[2]陈庆奎，王海峰，那丽春，等.图形处理器通用计算的研究综述[J]，黑龙江大学自然科学学报，2012，29（5）：672-679.

[3]陈茜.基于GPU的数字图像处理并行算法的研究[D].西安：中国科学院研究生院（西安光学精密机械研究所），2014.

[4]曾浩洋.图像处理的GPU加速技术的研究及实现[D].绵阳：西南科技大学，2019.

[5]张舒，褚艳利.GPU高性能运算之CUDA[M].北京：中国水利水电出版社，2009.

[6] Obukhov A.Kharlamov A.Discrete cosine transform for 8x8blocks with CUDA[R].Nvidia White Paper，2008.

[7]佘秋菊.基于DCT变换的JPEG图像压缩及其MATLAB实现[J].科技信息（学术研究），2008（36）：557-558.

[8]邹岩，杨志义，张凯龙.CUDA并行程序的内存访问优化技术研究[J].计算机测量与控制，2009，17（12）：2504-2506.

【通联编辑：谢媛媛】

作者简介：张金霜（1987-），男，广东茂名人，助教，硕士，主要研究方向为网络安全、智能优化算法;黄旭彬（1981-），广东高州人，讲师，硕士，主要研究方向为计算机应用技术。