第２７卷第８期 文章编号：１００６—９３４８（２０１０）０８—０２３１一ｏ４ 计算机仿真 ２０１０年８月 基于ＧＰＵ的小波变换 周 侃，阎文丽，甘斌，郝佳新 （中国人民６３９６１，北京１０００１２） 摘要：针对图形处理器研究问题，其中图形海量数据集的分析与处理，多用小波变换方法。但计算量大，难以适应实时性要 求。近年来图形处理器的性能大幅提高，其深度流水线和并行运算机制提高，为解决实时计算问题提供了良好的平台。在 研究小波变换矩阵形式及ＧＰＵ编程模型的基础上，提出了一种关于ＧＰＵ的小波变换方法，利用数组与纹理之间的对应关 系实现小波变换，将离散的数据点映射到纹理，将小波变换的计算影射为高维矩阵与向量间的乘积形式，并通过渲染到纹理 的形式取得中间结果。方法充分发挥了ＧＰＵ流水线的并行性优势，实验表明方法可有效减少计算时问，从而达到实时绘制 的要求。 关键词：图形处理器；小波变换；通用计算 中图分类号：ＴＰ３９１ 文献标识码：Ｂ Ｗａｖｅｌｅｔ Ａｎａｌｙｓｉｓ Ｂａｓｅｄ ｏｎ ＧＰＵ ＺＨＯＵ Ｋａｎ，ＹＡＮ Ｗｅｎ—ｌｉ，ＧＡＮ Ｂｉｎ，ＨＡ０ Ｊｉａ—ｘｉｎ （６３９６１ Ｕｎｉｔ，ＰＩ Ａ，Ｂｅｉｊｉｎｇ １０００１２，Ｃｈｉｎａ） ＡＢＳＴＲＡＣＴ：Ｗａｖｅｌｅｔ Ａｎａｌｙｓｉｓ ｈａｓ ｂｅｅｎ ｂｒｏｕｇｈｔ ｉｎｔｏ ｗｉｄｅ ｕｓｅ ｉｎ ｌａｒｇｅ ｓｃａｌｅ ｄａｔａ ｐｒｏｃｅｓｓｉｎｇ，ｂｕｔ ｉｔ ｃａｎ ｎｏｔ ｆｉｔ ｔｈｅ ｄｅｍａｎｄ ｏｆ ｒｅａｌ—ｔｉｍｅ ｓｉｍｕｌａｔｉｏｎ ｂｅｃａｕｓｅ ｔｈｅ ｑｕａｎｔｉｔｙ ｏｆ ｃｏｍｐｕｔｉｎｇ ｉｓ ｖａｓｔ．Ｉｎ ｒｅｃｅｎｔ ｙｅａｒｓ，ｔｈｅ ｃａｐａｂｉｌｉｔｙ ｏｆ ＧＰＵ ｉｓ ｉｍｐｒｏｖｅｄ ｇｒｅａｔｌｙ，ｔｈｅ ｐａｒａｌｌｅｌ ｐｉｐｅｌｉｎｅ ｐｒｏｃｅｓｓｏｒ ａｎｄ ｐｒｏｇｒａｍｍａｂｉｌｉｔｙ ｏｆ ＧＰＵ ｐｒｏｖｉｄｅ ａ ｓｏｌｕｔｉｏｎ ｔｏ ｔｈｉｓ ｑｕｅｓｔｉｏｎ．Ｗｅ ｐｒｅｓｅｎｔ ａ ｍｅｔｈｏｄ ｆｏｒ ｗａｖｅｌｅｔ ａｎａｌｙｓｉｓ ｂａｓｅｄ ｏｎ ＧＰＵ ａｎｄ ｍａｔｒｉｘ ａｒｉｔｈｍｅｔｉｃ ｏｆ ｗａｖｅｌｅｔ，ｗｈｉｃｈ ｕｓｅｓ ｔｈｅ ｍａｔｒｉｘ ａｎｄ ｔｅｘ— ｔｕｒｅ ｔｏ ｃｏｍｐｌｅｔｅ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎ，ｔｒａｎｓｆｏｒｍｓ ｔｈｅ ｄａｔａ ｔｏ ｔｅｘｔｕｒｅ，ｔｒａｎｓｆｏｒｍｓ ｔｈｅ ｗａｖｅｌｅｔ ａｎａｌｙｓｉｓ ｔｏ ｐｒｏｄｕｃｔ ｏｆ ｔｈｅ ｍａ— ｔｉｘ ａｎｄ ｖｅｃｔｒｏｒ，ａｎｄ ｕｓｅｓ ｔｈｅ ｒｅｎｄｉｎｇ ｔｏ ｇｅｔ ｔｈｅ ｒｅｓｕｌｔｓ ｏｆ ｃｏｍｐｕｔａｔｉｏｎ．Ｔｈｅ ｍｅｔｈｏｄ ｔａｋｅｓ ａｄｖａｎｔａｇｅ ｏｆ ＧＰＵ ｐｉｐｅｌｉｎｅ， ａｎｄ ｔｈｅ ｒｅｓｕｌｔｓ ｏｆ ｅｘｐｅｒｉｍｅｎｔ ｉｎｄｉｃａｔｅ ｔｈａｔ ｔｈｅ ｍｅｔｈｏｄ ｒｅｄｕｃｅｓ ｔｈｅ ｔｉｍｅ ｏｆ ｃａｌｃｕｌａｔｉｏｎ ｅｆｆｅｃｔｉｖｅｌｙ，ＳＯ ｔｈａｔ ｉｔ ａｃｈｉｅｖｅｓ ｔｈｅ ｄｅｍａｎｄ ｏｆ ｒｅａｌ～ｔｉｍｅ ｒｅｎｄｅｒ． ＫＥＹＷＯＲＤＳ：ＧＰＵ；Ｗａｖｅｌｅｔ ａｎａｌｙｓｉｓ；Ｇｅｎｅｒａｌ ｐｕｒｐｏｓｅ ｃｏｍｐｕｔａｔｉｏｎ １ 引言 小波理论是２Ｏ世纪８Ｏ年代后期发展起来的数学分支， 序力不从心，难以满足实时计算的要求。近年来图形处理器 （ＧＰＵ）的性能飞速发展，以ＧｅＦｏｒｅｅ ８８００系列显卡为例，其 拥有１２８个流处理器，核心频率达６１２ＭＨｚ，流处理器速度达 到了１．５０ＧＨｚ，可以每秒处理３３６亿个具有纹理的三角形 （Ｔｅｘｔｕｒｅｄ Ｔｒｉａｎｇｌｅ），而ＧＰＵ的可编程性也越来越好，这使得 近年来基于ＧＰＵ的通用计算（Ｇｅｎｅｒａｌ Ｐｕｒｐｏｓｅ ＧＰＵ，ＧＰＧＰＵ） 是目前国际公认的最新的空问（时间）一频率分析工具，由于 其“自适应性”和“数学显微镜性质”而成为许多学科共同关 注的焦点。由于它同时具有时域和频域的良好局部化性质， 所以可以随着信号不同频率成分在时空域采用疏密自动调 节的窗口分析信号特征，借助于小波分析，可以检测和提取 多源、多尺度、海量的空间数据集的基本特征，广泛应用于信 应运而生，它指的是利用ＧＰＵ来实现一般意义上的计算，而 不单纯是绘制¨ｌ 。ＧＰＵ用于通用计算的优势主要体现在 ｌ其高并行性的流水线结构上，通常情况下，ＧＰＵ拥有４条顶 点流水线和８条像素流水线可以对大量的计算进行并行处 理，再加上ＧＰＵ硬件内部实现了一些计算指令，所以ＧＰＵ在 众多计算工作中效率超过ＣＰＵ（如在向量计算方面效率为 ＣＰＵ的十倍）。 号处理…、图像处理 、三维地形分析 等方面。 但由于对大规模数据（例如：大规模ＤＥＭ数据）进行小 波变换的计算量庞大，这使得基于串行方式运行的ＣＰＵ程 基金项目：国家８６３高技术计划项目资助（２００６ＡＡ０１Ｚ３１９） 收稿日期：２００９—０５—１４修回日期：２００９—０７—０２ 本文在研究ＧＰＵ编程模型和小波变换矩阵形式的基础 上，提出了一种基于ＧＰＵ的小波变换方法，将离散小波变换 一２３１— 映射到纹理光栅化，充分发挥了ＧＰＵ的并行性，有效提高了 表示这种范围。 大规模数据小波变换的效率。 ６）计算的输出范围＝顶点坐标 光栅器产生片段，这些片段经过处理后变成输出像素。 ２ ＧＰＵ通用计算编程模型 输入顶点和顶点程序决定了应该生成哪些像素，因此顶点坐 在ＧＰＵ流水线中，顶点着色器（Ｖｅｒｔｅｘ Ｓｈａｄｅｒ）被用于控 标控制了计算的输出范围。 制渲染管线中的顶点几何变换，所以能够胜任各种几何计 算。而片段着色器（Ｆｒａｇｍｅｎｔ Ｓｈａｄｅｒ）被用于控制像素颜色 ３小波变换矩阵形式 和纹理映射等，所以具有较大容量的纹理空间，其对纹理数 由于ＧＰＵ是流处理器，需要将小波变换过程分割成一 一 １Ｊ Ｏ 据的操作可以广泛用于各种通用计算，如二维、三维、四维向 系列核心模块来处理作为纹理的数据流，而小波变换的矩阵 凡 量的加、减、乘、乘加、内积、最小值和最大值、取反、外积等， 形式可以天然地与纹理操作吻合，成为简单有效的核心模 ～ １ｊ Ｏ 它也是ＧＰＵ用于通用计算的主要途径。ＧＰＵ通用计算的实 块。利用多分辨率分析（ＭＲＡ）理论，可以构造离散小波基 ． ．ｒＬ 质其实是绘制的过程，以下介绍ＧＰＵ通用计算的编程模型： 计算框架，并计算出不同小波基的滤波器系数，以ｂｉ一 一 Ｏｊｏｒ ３．３小 １）ＣＰＵ数组＝ＧＰＵ纹理 波为例，可以得到其低通滤波器系数和高通滤波器系数： ．ｒｌ ｒＬ 任何ＣＰＵ上使用的数组都可以用ＧＰＵ上的纹理来代 ｛ｈ ｝＝｛０．０４６９，一０．１４０６，一０．１０９４，０．７０３１，０．７０３１，一 ７ １ ］●Ｊ １ｊ 替，并通过片段着色器对其进行处理。 一０．１０９４，一０．１４０６，０．０４６９｝ ２）ＣＰＵ内循环＝ＧＰＵ片段程序 ｛ｇ ｝＝｛０，０，０．６３６９，一１．９１０７，１．９１０７，一０．６３６９，０，一 ０ ¨ 在ＣＰＵ上用循环语句（如ｆ０ｒ循环）来操作数组中的元 ０｝ 素，而在ＧＰＵ片段程序中也可以编写相似的指令来处理 由离散小波分解公式可知： ～Ｏ １Ｊ 纹理。 ７ ［ｋ］＝∑ｈ［ｎ］ ｎ＋２ ］ 一 Ｏ ３）反馈：渲染到纹理 （１） ｎ ０ 在ＣＰＵ上得到计算结果是方便的，因为它的存储器可 这里将低通滤波器的下标写到中括号里。式（１）是一个 一 ０ 以在程序中的任何地方进行读写，而ＧＰＵ为了实现计算反 迭代过程，给定一个初始序列｛ｃｊ｝，就可以计算出｛ｃｊ一 ｝， 馈，必须使用渲染到纹理把片段程序的结果写入存储器。 ｛ｃｊ一 ｝，…，由于通常的迭代是递增过程，所以将式（１）写成 ４）计算的调用＝几何体光栅化 ７ ＧＰＵ通用计算的实质是绘制的过程，所以为了调用计 如下形式：ｃ，＝＋Ｉ［ ］＝∑ ［ｎ］ｎ ０ ［ｎ＋２｜ｊ｝］ 算，只需要绘制几何体。 假设对输入信号 ｔ）进行整数点采样，得到序列：｛．．・ ５）计算的输入范围＝纹理坐标 ［０］ １］ ２］，… ｎ一１］，…｝，令初始迭代序列ｃ。［ｋ］：， 任何计算都有一个输入范围。ＧＰＵ以纹理坐标的形式 ［ｋ］，（ｋ Ｚ），于是小波变换公式可以写成如下的矩阵形式： －●● ｃｊ＋ ［０］ ｃｉ［Ｏ］ ｃｊ＋ ［１］ ｃｊ［１］ （２） ●●● ｃｊ［ｎ一１］ 令低通变换矩阵Ｐ为： Ｐ＝ （３） ＝（…，ｃｊ［Ｏ］，ｃｊ［１］，…，ｃｉ［ｎ一１］，…） 于是Ｍａｌｌａｔ低通分解公式可写成如下矩阵形式： ＋ ＝ＰｃＪ，（ ＝０，１，…，Ｌ一１） 表示小波分解的次数，同理，得到高通变换矩阵为： —－－——２３２－－－—— …ｇ［０］ｇ［１］ ０ ０ … ｇ［７］０］ｇ［１］ ０ ０ … … ０ ・ （４） …Ｑ＝ … ｇ［７］０ … ｇ［ｏ］ｇ［１］…ｇＥ７］・ ＝（…，Ｃｏｌ， ［１］，…， ［ｎ一１］，…） 重构出原始信号的每一行。这个过程正好是二维小波分解 矩阵形式的逆过程。 ｒ Ｃｏ＝ｆ ｛ｑ＝ｐ ．．，（ ＝１，２，…， ） （５） ＝Ｑ 一・ 二维小波分解矩阵形式如图１所示：提取原始数据的每 一行分别与高通滤波器系数矩阵、低通滤波器系数矩阵相 乘，得到结果后之后再提取每一列继续和两个滤波器系数矩 阵相乘，最后得到小波分解的四个低频系数矩阵。 图１二维小波分解矩阵形式 二维小波重构也可以得到相应的矩阵形式。在一维小 波变换中，由双尺度方程可以及ｇ ＝（一１） ｈ 一 的关系可 以得到： ∑ｈｋ＝０ ：０ … … 一 ０ 尸ｒ尸＋Ｑ Ｑ＝ ＝Ｅ ｋ＝０ Ｅ ｈ：０ … … … ０ （６） 其中Ｅ为单位阵，在上式两边同乘信号Ｃ ，得到： 一。＝（Ｐ Ｐ＋Ｑ Ｑ） 一。＝ｐｒＰ 一 ＋Ｑ Ｑ 一。 再由式（５）可知： 一１＝Ｐ Ｇ，＋Ｑ Ｄ， （７） 利用式（７）可以得到二维小波重构的矩阵方法（如图２ 所示）：首先分别提取低频分量和三个高频分量矩阵的每一 列，分别与相应的矩阵（Ｐ 和Ｑ ）相乘，之后相加，在得到的 结果中提取每一行，分别与相应的矩阵相乘之后相加，便可 图２二维小波重构矩阵形式 ４基于ＧＰＵ的小波变换实现 利用ＧＰＵ处理海量数据的小波变换的好处是很明显 的，由式（２）可以看出，以ｎ×ｎ的二维离散小波变换分解为 例，进行一次小波变换相当于利用高通滤波器和低通滤波器 对数据块进行行方向与列方向一共４次滤波，相当于２ｎ 次 两个ｎ维向量的张量积，也相当于大约２ｎ 次乘法和２ｎ 次 加法，这些高密集度的计算在ＣＰＵ上是利用顺序执行的循 环语句运行的，而对于流水线结构的ＧＰＵ而言，将并行处理 作为向量的纹理中的每个纹素，也就是说在一个时钟周期处 理多个纹理元素，因此大大加快了计算速度。另外，显卡内 存接口为２５６位（以ＧｅＦｏｒｃｅＦＸ为例），大于ＣＰＵ的３２位接 口，所以计算带宽要明显大于ＣＰＵ。最后当绘制时，由于计 算数据已在显存中，减少了内存与显存的传输（这种传输相 对较慢），节省了ＣＰＵ资源，从而可以使ＣＰＵ有空做更多的 控制操作。 基于ＧＰＵ的小波变换方式，利用小波变换的矩阵化方 法将小波变换的计算影射为一个高维矩阵与向量间的乘积 形式，而中间结果为一些向量，利用数组与纹理之间的的对 应关系实现小波变换。小波分解时的基本思想如下： 首先将式（５）中的向量影射到纹理，自然想到ＧＰＵ的一 维纹理，但由于ＧＰＵ上的一维纹理长度是有限的，最长为 ４０９６，所以最多只能实现一个４０９６长度的向量，另外，由于 作为计算结果的向量需要被渲染，而ＧＰＵ在生成相同片段 个数的前提下，渲染二维纹理比渲染一维纹理快的多。所 以，可以将一维向量映射到二维纹理上进行处理，为了进一 步降低内部表示数据的大小，将连续的４个向量中的元素压 缩为一个ＲＧＢＡ纹理元素，如图３所示。 ．．．——２３３．．．—— 兰！Ｉ竺 ｌ！ ｌ竺 ｌ！！ｌ！ ｌ！ ｌ！！ｌ竺！ｌ兰！ｌ！ ！ ！ ｆ ｑ Ｒ Ｇ Ｂ Ａ 加，对得到的两个矩阵继续提取每一行，再次计算矩阵和Ｐ 、 Ｒ Ｇ Ｂ Ａ Ｑ 的乘法并将结果相加，最后得到重构出的原始信号。 鹾盎 图３将向量表示为二维纹理 ５ 实验结果 利用小波变换的ＣＰＵ方法和ＧＰＵ方法对不同大小的数 据块进行分解和重构（实验平台：ＣＰＵ：Ｉｎｔｅｌ Ｃｏｒｅ ２ Ｄｕｏ Ｅ４５００ ２．２ＧＨ，内存：１ＧＢ ＤＤＲ２ ８００，显卡：８６００ＧＴ ５１２ＭＢ显 存），两种方法所用时问如图４所示。 而对于式（５）中的滤波器系数矩阵，则可以先将矩阵拆 成ｎ个向量，对于每个向量再利用上边的二维纹理的形式来 存储。 对于式（２）的矩阵与向量之间乘积运算，可以先将其转 化为ｎ次两个向量的乘积，而对于作为纹理的向量间的乘 积，则利用绘制双重纹理，首先把渲染区域设置为与表示向 量的二维纹理同样的屏幕尺寸，然后将纹理赋给一个覆盖整 个屏幕（渲染区域）的方块，顶点到达光栅化阶段时，光栅器 为每个向量元素（纹理）生成一个片段，每个片段在经过片段 程序时，可将两个向量元素取回并作乘法运算，然后将结果 写回到渲染目标。下面以ｎ×ｎ二维离散小波分解为例，将 该方法分为如下几个步骤： １）利用事先求得的对称小波滤波器组系数｛ｈ ｝和｛ｇ ｝ 构造出相应的两个滤波器系数矩阵Ｐ和Ｑ（如式（３）、（４）的 形式），对于矩阵中的每一行，都将其作为一张二维纹理进行 保存，另外将待变换的数据（ｎ×ｎ）的每一行每一列分别保存 为二维纹理。 ２）创建一张新的纹理用来保存计算结果。 ３）设置绘制区域大小为ｎ×ｎ个像素，并且绘制一个１： １全屏幕矩形。 ４）对二维数据进行行方向上的低通滤波，即将Ｐ的第一 行所对应的二维纹理和二维数据的第一行所对应的二维纹 理作为双重纹理赋予矩形，令双重纹理计算方式为乘法。利 用渲染到纹理取回渲染结果，对于每次得到的结果纹理的各 像素灰度值相加，得到最终结果。 ５）对于Ｐ的每一行，执行步骤４的计算，将得到的ｎ／２ 个最终结果，将这些结果按顺序排列，便得到了二维数据第 一行的低通滤波下采样结果。 ６）同理，对二维数据的每一行进行步骤４）、步骤５）的计 算，得到二维数据经过行方向上低通滤波下采样结果。 ７）利用与步骤３）、４）、５）、６）相同的方式，将Ｑ与二维数 据进行双重纹理计算，可以得到二维数据经过行方向上高通 滤波下采样结果。 ８）将得到的行方向上低通和高通滤波下采样结果按列 划分为二维纹理，利用低通滤波器和高通滤波器对其进行滤 波，利用步骤１）～７），最终得到数据块小波变换后的４幅系 数子图。 重构与分解的原理是一致的，首先将已知的低频小波系 数矩阵和三个高频小波系数矩阵进行按列的提取，并且利用 双重纹理计算这四个矩阵和Ｐ。。’、Ｑ 的乘法，之后将结果相 ．．．——２３４．．．—— ｇ 莒 图４ ＣＰＵ和ＧＰＵ小波变换运算时间比较 从实验数据可以看出，ＧＰＵ的计算时间要远小于ＣＰＵ 的计算时间，并且ＣＰＵ的运算时间随着数据量的增长急剧 增加，而ＧＰＵ的耗时增长较少，从耗时比可以看出，随着计 算量的增大，ＧＰＵ表现出更优越的性能。 参考文献： ［１］张仁辉，杜民．小波分析在信号去噪中的应用［Ｊ］．计算机仿 真，２００５，２２（８）：６９—７２． ［２］李云，刘学诚．小波变换在图像处理中的应用［Ｊ］．计算机仿 真，２００８，２５（６）：１９５—１９７． ［３］ 常占强，吴立新．用离散小波变换对格网ＤＥＭ数据压缩中的 关键技术［Ｊ］．系统仿真学报，２００８，２０（１５）：３９５５—３９６２． ［４］ Ｗｕ Ｅｎｈｕａ，Ｌｉｕ Ｙｏｕｑｕａｎ．Ｅｍｅ ｎｇ ｔｅｃｈｎｏｌｏｇｙ ａｂｏｕｔ ＧＰＧＰＵ 『Ｃ］．Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ｔｈｅ ２００８ ＩＥＥＥ Ａｓｉａ Ｐａｃｉｉｆｃ Ｃｏｎｆｅｒｅｎｃｅ．Ｍａ— ｃａｏ。Ｃｈｉｎａ．２００８．６１８—６２２． ［５］Ｊ Ｏｗｅｎｓ，Ｄ Ｌｕｅｂｋｅ，Ｎ Ｇｏｖｉｎｄａｒａｊｕ，Ｍ Ｈａｒｉｒｓ．Ａ ｓｕｒｖｅｙ ｏｆｇｅｎｅｒ－ ａｌ—ｐｕｒｐｏｓｅ ｃｏｍｐｕｔａｔｉｏｎ ｏｎ ｇｒａｐｈｉｃｓ ｈａｒｄｗａｒｅ［Ｊ］．Ｃｏｍｐｕｔｅｒ Ｇｒａｐｈｉｃｓ Ｆｏｒｕｍ，２００７，２６（１）：８Ｏ一１ １３． ［作者简介］ 周侃（１９８３一），男（汉族），吉林春市人，硕 士研究生，助理工程师，主要研究领域为多媒体与 虚拟现实。 阎文丽（１９７２一），女（汉族），山西省清徐县人，硕 士研究生，高级工程师，主要研究领域为系统仿真。 甘斌（１９７９一），男（汉族），江西九江市人，硕士研究生，工程师， 主要研究领域为系统仿真。 郝佳新（１９８０一），男（汉族），河北三河市人，硕士研究生，工程师， 主要研究领域为系统仿真。
