基于CUDA技术的离散小波变换算法研究与实现

摘  要：针对离散小波变换过程比较耗时、不利于实际工程应用的问题，提出利用基于GPU平台的CUDA技术对小波变换算法做并行化改造，从而提高算法执行效率。该文分析了小波Mallat算法并行化的可行性，并详细介绍了算法的改造过程。实验表明，基于GPU/CUDA技术的并行小波Mallat算法，相较于串行小波变换算法，执行速度最高提升了50余倍，且算法效率与计算量成正向关系。

关键词：CUDA;并行程序设计;离散小波变换;图像压缩

中图分类号：TN911;TP399       文献标识码：A 文章编号：2096-4706（2020）17-0072-05

Abstract：In view of the time-consuming process of discrete wavelet transform and not conducive to practical engineering applications，it is proposed to use GPU-based CUDA technology to parallelize the wavelet transform algorithm，thereby improving the efficiency of algorithm execution. This article analyzes the feasibility of the parallelization of the wavelet Mallat algorithm，and introduces the algorithm transformation process in detail. Experiments show that the parallel wavelet Mallat algorithm based on GPU/CUDA technology，compared with the serial wavelet transform algorithm，the execution speed is up to 50 times faster，and the algorithm efficiency has a positive relationship with the amount of calculation.

Keywords：CUDA;parallel programming;discrete wavelet transform;image compression

0  引  言

小波分析廣泛应用于信号分析、图像识别、计算机视觉和数据压缩等，JPEG2000标准的出现则标志着小波变换在图像压缩方面走向成熟，Mallat算法将小波变换真正应用于实际，该算法采用滤波器执行离散小波变换，但对于数据量较大的图像信息，小波变换的过程耗时较长。为了提升小波变换算法的执行效率，可考虑将传统CPU上的串行结构算法改成基于CPU/GPU异构算法，通过GPU内部大量独立的计算单元，实现计算多线程并发执行[1-3]。

为方便作者在后续的图像分析研究工作中，使用高效的离散小波变换技术，及在自研的工程项目中便捷地调用小波变换算法，因此需要自己重新编写小波变换程序。本文将结合小波Mallat算法的数学理论与程序代码逻辑，论证该算法具备并行化改造的可行性;再从程序设计优化的角度，利用CUDA技术改写小波Mallat算法中的耗时部分，通过利用GPU强大的计算资源来提升算法效率。

1  CUDA技术

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

在CUDA架构下，程序被分成两个部分：host端程序和device端程序。host端程序运行在CPU上，device端程序又叫kernel，运行在GPU上。通常先由host端程序准备数据，对设备进行必要的初始化，把数据从内存传送到显存中，执行device端程序，待GPU运行完后，host端程序再从显存中把结果取回内存中。

CUDA的线程模型分三个层次：线程（thread）、线程块（block）、线程网格（grid）。其中，kernel以grid的形式组织，每个线程网格由若干个block组成，而每个线程块又由若干个thread组成。block内的thread之间是并行执行的，而grid内的block之间也是并行执行的，grid与grid之间是串行执行的，即GPU同时只能处理一个grid[4]。

2  离散小波变换算法

离散小波及离散小波变换的表达式为[5]：

其中a为尺度因子，b为平移因子，且a>1，b∈R;m，n∈Z2，m为频率范围指数，n为时间步长变化指数;ψm，n（t）为基本小波或母小波;f（t）为原函数，〈〉为平方可积函数空间L2（R），而〈f，ψm，n〉为小波变换。

离散小波变换是一种时频分析，它从集中在某个区间上基本函数开始，以规定步长向左或者向右移动基本波形，并用尺度因子a扩张或收缩来构造函数的系数。

如果尺度以二进方式离散变化，可得小波函数为：

其中参数j决定伸缩，而k确定平移幅度，且j，k∈Z。

由于图像信息一般是二维的，因此需要将小波理论从一维推广至二维，二维基本小波如式（4）所示。

假设一幅M×M像素的图像f1（x，y），根据Mallat算法，对于j=0，则源图像的尺度为20=1。j值每增大1都使尺度加倍，分辨率减半。每次隔行隔列采样，图像被分解为4个大小为原来尺寸四分之一的子频带区域，这4个子频带区域的每一个都是由源图像与一个小波基函数内积后，再经过x和y方向都进行2倍的间隔抽样而产生的。因此，对于第一个层次j=1，可写成：

其中f20（m，n）为频带区域，其下标为尺度，上标为索引而非指数。

对于后继的层次（j>1），（x，y）都以完全相同的方式分解而构成4个尺度上的更小的图像。由此可知，执行小波变换的有效方法是使用滤波器，将内积写成卷积的形式：

因为尺度函数和小波函数都是可分离的，所以每个卷积都可以分解成在 （x，y）的行和列上的一维卷积。由式

（6）可知，实现Mallat小波分解的核心是将待处理数据与滤波器进行卷积，最后得到一系列小波系数[6]。

3  基于CUDA实现Mallat算法加速

3.1  算法整体设计

要利用GPU来加速变换过程，则必须结合GPU的特性来优化算法的效率。算法实现的总体流程如图1所示。

算法从host端开始，首先在内存中分配两个具有与被处理图像像素数据相同大小的空间，分别用于存放被处理图像的源数据和处理结果，并把图像像素数据读入到内存。图像要在GPU中进行处理，因此需要把图像像素数据复制到显存中。在显存中申请两个具有与被处理图像像素数据相同大小的空间，分别用于存放图像源数据和处理结果。其中，存放图像源数据的空间采用纹理存储器，存放处理结果的空间采用全局存储器。

每一个线程处理一个像素点的卷积运算。小波变换完成后，将结果写入到全局存储器中，同样是每个线程写入一个系数，至此device端的工作完成。最后，host端进行数据量化和小波重構，并显示最终结果，整个算法完成。

3.2  源程序分析

本文采用的小波Mallat算法的伪代码为：

//卷积操作

void Convolution {

for（..;..;..） {

for（..;..;..） {

行（列）前半部分时域数据 += 时域数据*分解尺度函数;

行（列）后半部分时域数据 += 时域数据*分解母函数;

}

}

}

//离散小波变换过程

void Wavelet（） {

for（n=0;n

Height=lHeight>>2;     //图像低频信息的高度

Width=lWidth>>2;       //图像低频信息的宽度

for（i = 0; i < Height; i++） {

for（j=0 ; j< Width ; j++）

{ 获取行数据 }

Convolution（…）;   // 对行方向进行卷积运算

}

for（i = 0; i < Width; i++） {

for（j=0;j< Height;j++）

{ 获取列数据 }

Convolution（…）;   // 对列方向进行卷积运算

}

} //小波变换结束

}

循环控制小波变换的层数，每循环一次，则进行一次完整的行卷积运算和列卷积运算，同时小波变换处理的范围缩小为原来的1/4。此时图像的低频有效信息集中在图像的左上角，而其他区域为高频细节信息，不属于下一次循环小波处理的区域。

由上述代码可知，该算法耗时的部分出现在图像的遍历和卷积过程，改进的关键在于如何将该部分串行代码变成运行在GPU上并行代码。

并行算法设计思路为：

（1）串行算法中，通过两重循环，实现图像的遍历。针对这个操作，在CUDA中，可以将源图像数据映射到纹理存储器（这过程称为纹理映射），利用纹理存储器支持二维寻址的特性，实现对像元的快速访问。

（2）串行算法中，每个像素点要完成卷积运算，即使每个像素点的计算相对独立，也必须依次进行。针对这个操作，在CUDA中，可以开辟多个线程，每个线程处理一个像素点，独立完成卷积。

3.3  并行算法设计

3.3.1  线程划分

在GPU的SM（Streaming Multiprocessor，多处理器）中，线程被组织成warp（线程束）执行，一个warp中包含32个线程，所以一个block中的线程数最好是32的倍数，这里选择一个block中申请16×16个thread，即定义dim3 threadBlock（16，16）。把整个图像作为一个grid，整个grid包含和图像像素点相等的线程数，则block的数量为（width/threadBlock.x）×（height/threadBlock.y），其中，width和height分别表示图像的宽和高。线程划分如图2所示。

3.3.2  存储器分配

为了存放源图像数据和处理后的图像数据，必须开辟两个和源图像数据相同大小的存储空间。

首先，存放源图像数据可采用CUDA数组，这可方便地将数据映射到纹理存储器。纹理数据按16的倍数对齐，符合“合并访问”的要求，能提高访存效率。由于在做列卷积的过程中，需要用到行卷积后的数据，为此要保存原来的行卷积结果，并更新纹理存储器中数据。纹理存储器是个只读存储器，不支持数据写入，但可以通过更新与之绑定的CUDA数组中的值，来达到更行纹理存储器中的数据的目的。这个过程中，因为数据仍然在device端进行传递，耗时不长。

其次，存放小波变换处理后的数据，可采用全局存储器。全局存储器使用的是普通显存，整个grid的所有线程都能访问，由于没有片内缓存，访问速度慢，但通过合并访问，可以隐藏延时。

此外，卷积过程中用到的“分解尺度函数”和“分解母函数”初始化后不再发生改变，故可将它们存放在常数存储器中。常数存储器大小为64 KB，是只读的地址空间，其数据位于显存，但拥有缓存加速。在CUDA程序中常数存储器一般用于存储需要频繁访问的只读参数[7]。

3.3.3  并行算法实现

host端控制小波变换的层数，device端执行一次小波變换。host端关键伪代码为：

for（int n=0;n

width=lWidth>>n    ;  //通过移位操作，控制小波变化的范围

height=lHeight>>n ;  //通过移位操作，控制小波变化的范围

//根据小波变换的层次，不断调整block的数量

dim3 blockGrid（width/threadBlock.x， height/threadBlock.y）;

convolutionRow<<>>（.. ..）; //行方向卷积

cudaMemcpyToArray（cu_array，0，0，d_Result，size，cudaMemcpyDeviceToDevice）; //更新纹理数据

convolutionCol<<>>（.. ..）; //列方向卷积

cudaMemcpyToArray（cu_array，0，0，d_Result，size，cudaMemcpyDeviceToDevice）; //更新纹理数据

}

因为每次小波变换处理的区域都在发生改变，为了在device端更方便地控制线程操作，故每次循环都重新定义grid，修改grid中block的数量。

接下来，调用device端的行（列）卷积运算，进行小波变化，每完成一次行（列）卷积，都必须更新纹理数据，以供下一次循环数据调用。device端小波变换伪代码为：

//行方向上卷积

if（x

for（int j=0; j< 某参数; j++）{

m=（某计算）% width;  // width是小波变换作用的区域

tmp += tex2D（texData， （（float）m + 0.5f） ， iy） *LF[j]; //时域与分解尺度函数相乘

}

f[y*lwidth+x] = tmp; // lwidth是原图像的宽度，f存放结果

}else if（x

for（int j=0; j< 某参数; j++）{

m=（某计算）% Width; // width是小波变换作用的区域

tmp += tex2D（texData， （（float）m + 0.5f） ， iy） *HF[j]; //时域与分解母函数相乘

}

f[y*lwidth+x] = tmp; // lwidth是原图像的宽度，f存放结果

}

上述代码的执行过程是并行的，即每个线程会并发地去执行相应的计算，然后将计算的结果保存在全局存储器中。列方向卷积同理，不再列举。

由上述分析可见，算法的时间复杂度由原来的四重循环变为两重循环，而两重循环的基数非常小，所以大大地缩短了算法运行时间。

4  实验结果分析

本程序目的在于对比小波变换在CPU和GPU上的运行效率，通过调用同一幅图片进行小波变换，取得各自的执行时间，最后算出加速比。

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

实验程序采用Daubechies D4小波，变换层次为3。所有数据为5次计算时间平均值，CUDA计算时间包括了数据在host端（主机）与device端（显卡存储器）中交换的时间。

经过测试，小波变换在CPU和GPU上执行的效率对比如表1所示。从表中数据可以看出，应用CUDA技术加速小波变换，均达到两个数量级的加速效果。

由图3可知，随着图像大小的增加，CPU上小波变换所需的时间急剧上升，而GPU上执行所需时间则增长缓慢。由图4可见，图像大小与加速比总体呈正向关系。

5  结  论

本文重点分析了小波Mallat算法并行化的可行性，并成功利用CUDA技术实现了高速并行小波变换算法。实验表明，相较于CPU上的串行小波变换算法，本文算法的执行效率提升了数十倍，最高达到了约51.6倍，且算法效率与计算量成正向关系。通过合理设计CUDA程序，将耗时的串行计算过程并行化，充分利用GPU强大的计算资源，能最大限度地缩短执行时间，使得小波技术更好地应用于实际工程系统中。

参考文献：

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

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

[3] 白兆峰.基于GPU的JPEG2000图像压缩编码技术的研究 [D].西安：中国科学院大学（中国科学院西安光学精密机械研究所），2016.

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

[5] 陈天华.数字图像处理 [M].北京：清华大学出版社，2007：132-133.

[6] 韩志伟，刘志刚，鲁晓帆，等.基于CUDA的高速并行小波算法及其在电力系统谐波分析中的应用 [J].电力自动化设备，2010，30（1）：98-101+105.

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

作者简介：张金霜（1987—），男，汉族，广东茂名人，教师，毕业于华南农业大学，硕士研究生，研究方向：计算机教育。