基于CUDA技術(shù)的離散小波變換算法研究與實現(xiàn)

摘? 要：針對離散小波變換過程比較耗時、不利于實際工程應(yīng)用的問題，提出利用基于GPU平臺的CUDA技術(shù)對小波變換算法做并行化改造，從而提高算法執(zhí)行效率。該文分析了小波Mallat算法并行化的可行性，并詳細介紹了算法的改造過程。實驗表明，基于GPU/CUDA技術(shù)的并行小波Mallat算法，相較于串行小波變換算法，執(zhí)行速度最高提升了50余倍，且算法效率與計算量成正向關(guān)系。

關(guān)鍵詞：CUDA;并行程序設(shè)計;離散小波變換;圖像壓縮

中圖分類號：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? 引? 言

小波分析廣泛應(yīng)用于信號分析、圖像識別、計算機視覺和數(shù)據(jù)壓縮等，JPEG2000標準的出現(xiàn)則標志著小波變換在圖像壓縮方面走向成熟，Mallat算法將小波變換真正應(yīng)用于實際，該算法采用濾波器執(zhí)行離散小波變換，但對于數(shù)據(jù)量較大的圖像信息，小波變換的過程耗時較長。為了提升小波變換算法的執(zhí)行效率，可考慮將傳統(tǒng)CPU上的串行結(jié)構(gòu)算法改成基于CPU/GPU異構(gòu)算法，通過GPU內(nèi)部大量獨立的計算單元，實現(xiàn)計算多線程并發(fā)執(zhí)行[1-3]。

為方便作者在后續(xù)的圖像分析研究工作中，使用高效的離散小波變換技術(shù)，及在自研的工程項目中便捷地調(diào)用小波變換算法，因此需要自己重新編寫小波變換程序。本文將結(jié)合小波Mallat算法的數(shù)學理論與程序代碼邏輯，論證該算法具備并行化改造的可行性;再從程序設(shè)計優(yōu)化的角度，利用CUDA技術(shù)改寫小波Mallat算法中的耗時部分，通過利用GPU強大的計算資源來提升算法效率。

1? CUDA技術(shù)

CUDA是由顯卡廠商NVIDIA推出的通用并行計算框架，允許設(shè)計人員使用C/C++語言和CUDA擴展庫的形式編寫程序，避免了傳統(tǒng)方法中將通用計算轉(zhuǎn)換為圖形處理的過程，大大降低程序設(shè)計的難度和系統(tǒng)維護成本。

在CUDA架構(gòu)下，程序被分成兩個部分：host端程序和device端程序。host端程序運行在CPU上，device端程序又叫kernel，運行在GPU上。通常先由host端程序準備數(shù)據(jù)，對設(shè)備進行必要的初始化，把數(shù)據(jù)從內(nèi)存?zhèn)魉偷斤@存中，執(zhí)行device端程序，待GPU運行完后，host端程序再從顯存中把結(jié)果取回內(nèi)存中。

CUDA的線程模型分三個層次：線程（thread）、線程塊（block）、線程網(wǎng)格（grid）。其中，kernel以grid的形式組織，每個線程網(wǎng)格由若干個block組成，而每個線程塊又由若干個thread組成。block內(nèi)的thread之間是并行執(zhí)行的，而grid內(nèi)的block之間也是并行執(zhí)行的，grid與grid之間是串行執(zhí)行的，即GPU同時只能處理一個grid[4]。

2? 離散小波變換算法

離散小波及離散小波變換的表達式為[5]：

其中a為尺度因子，b為平移因子，且a>1，b∈R;m，n∈Z2，m為頻率范圍指數(shù)，n為時間步長變化指數(shù);ψm，n（t）為基本小波或母小波;f（t）為原函數(shù)，〈〉為平方可積函數(shù)空間L2（R），而〈f，ψm，n〉為小波變換。

離散小波變換是一種時頻分析，它從集中在某個區(qū)間上基本函數(shù)開始，以規(guī)定步長向左或者向右移動基本波形，并用尺度因子a擴張或收縮來構(gòu)造函數(shù)的系數(shù)。

如果尺度以二進方式離散變化，可得小波函數(shù)為：

其中參數(shù)j決定伸縮，而k確定平移幅度，且j，k∈Z。

由于圖像信息一般是二維的，因此需要將小波理論從一維推廣至二維，二維基本小波如式（4）所示。

假設(shè)一幅M×M像素的圖像f1（x，y），根據(jù)Mallat算法，對于j=0，則源圖像的尺度為20=1。j值每增大1都使尺度加倍，分辨率減半。每次隔行隔列采樣，圖像被分解為4個大小為原來尺寸四分之一的子頻帶區(qū)域，這4個子頻帶區(qū)域的每一個都是由源圖像與一個小波基函數(shù)內(nèi)積后，再經(jīng)過x和y方向都進行2倍的間隔抽樣而產(chǎn)生的。因此，對于第一個層次j=1，可寫成：

其中f20（m，n）為頻帶區(qū)域，其下標為尺度，上標為索引而非指數(shù)。

對于后繼的層次（j>1），（x，y）都以完全相同的方式分解而構(gòu)成4個尺度上的更小的圖像。由此可知，執(zhí)行小波變換的有效方法是使用濾波器，將內(nèi)積寫成卷積的形式：

因為尺度函數(shù)和小波函數(shù)都是可分離的，所以每個卷積都可以分解成在 （x，y）的行和列上的一維卷積。由式

（6）可知，實現(xiàn)Mallat小波分解的核心是將待處理數(shù)據(jù)與濾波器進行卷積，最后得到一系列小波系數(shù)[6]。

3? 基于CUDA實現(xiàn)Mallat算法加速

3.1? 算法整體設(shè)計

要利用GPU來加速變換過程，則必須結(jié)合GPU的特性來優(yōu)化算法的效率。算法實現(xiàn)的總體流程如圖1所示。

算法從host端開始，首先在內(nèi)存中分配兩個具有與被處理圖像像素數(shù)據(jù)相同大小的空間，分別用于存放被處理圖像的源數(shù)據(jù)和處理結(jié)果，并把圖像像素數(shù)據(jù)讀入到內(nèi)存。圖像要在GPU中進行處理，因此需要把圖像像素數(shù)據(jù)復(fù)制到顯存中。在顯存中申請兩個具有與被處理圖像像素數(shù)據(jù)相同大小的空間，分別用于存放圖像源數(shù)據(jù)和處理結(jié)果。其中，存放圖像源數(shù)據(jù)的空間采用紋理存儲器，存放處理結(jié)果的空間采用全局存儲器。

每一個線程處理一個像素點的卷積運算。小波變換完成后，將結(jié)果寫入到全局存儲器中，同樣是每個線程寫入一個系數(shù)，至此device端的工作完成。最后，host端進行數(shù)據(jù)量化和小波重構(gòu)，并顯示最終結(jié)果，整個算法完成。

3.2? 源程序分析

本文采用的小波Mallat算法的偽代碼為：

//卷積操作

void Convolution {

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

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

行（列）前半部分時域數(shù)據(jù) += 時域數(shù)據(jù)*分解尺度函數(shù);

行（列）后半部分時域數(shù)據(jù) += 時域數(shù)據(jù)*分解母函數(shù);

}

}

}

//離散小波變換過程

void Wavelet（） {

for（n=0;n

Height=lHeight>>2;? ? ?//圖像低頻信息的高度

Width=lWidth>>2;? ? ? ?//圖像低頻信息的寬度

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

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

{ 獲取行數(shù)據(jù) }

Convolution（…）;? ?// 對行方向進行卷積運算

}

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

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

{ 獲取列數(shù)據(jù) }

Convolution（…）;? ?// 對列方向進行卷積運算

}

} //小波變換結(jié)束

}

循環(huán)控制小波變換的層數(shù)，每循環(huán)一次，則進行一次完整的行卷積運算和列卷積運算，同時小波變換處理的范圍縮小為原來的1/4。此時圖像的低頻有效信息集中在圖像的左上角，而其他區(qū)域為高頻細節(jié)信息，不屬于下一次循環(huán)小波處理的區(qū)域。

由上述代碼可知，該算法耗時的部分出現(xiàn)在圖像的遍歷和卷積過程，改進的關(guān)鍵在于如何將該部分串行代碼變成運行在GPU上并行代碼。

并行算法設(shè)計思路為：

（1）串行算法中，通過兩重循環(huán)，實現(xiàn)圖像的遍歷。針對這個操作，在CUDA中，可以將源圖像數(shù)據(jù)映射到紋理存儲器（這過程稱為紋理映射），利用紋理存儲器支持二維尋址的特性，實現(xiàn)對像元的快速訪問。

（2）串行算法中，每個像素點要完成卷積運算，即使每個像素點的計算相對獨立，也必須依次進行。針對這個操作，在CUDA中，可以開辟多個線程，每個線程處理一個像素點，獨立完成卷積。

3.3? 并行算法設(shè)計

3.3.1? 線程劃分

在GPU的SM（Streaming Multiprocessor，多處理器）中，線程被組織成warp（線程束）執(zhí)行，一個warp中包含32個線程，所以一個block中的線程數(shù)最好是32的倍數(shù)，這里選擇一個block中申請16×16個thread，即定義dim3 threadBlock（16，16）。把整個圖像作為一個grid，整個grid包含和圖像像素點相等的線程數(shù)，則block的數(shù)量為（width/threadBlock.x）×（height/threadBlock.y），其中，width和height分別表示圖像的寬和高。線程劃分如圖2所示。

3.3.2? 存儲器分配

為了存放源圖像數(shù)據(jù)和處理后的圖像數(shù)據(jù)，必須開辟兩個和源圖像數(shù)據(jù)相同大小的存儲空間。

首先，存放源圖像數(shù)據(jù)可采用CUDA數(shù)組，這可方便地將數(shù)據(jù)映射到紋理存儲器。紋理數(shù)據(jù)按16的倍數(shù)對齊，符合“合并訪問”的要求，能提高訪存效率。由于在做列卷積的過程中，需要用到行卷積后的數(shù)據(jù)，為此要保存原來的行卷積結(jié)果，并更新紋理存儲器中數(shù)據(jù)。紋理存儲器是個只讀存儲器，不支持數(shù)據(jù)寫入，但可以通過更新與之綁定的CUDA數(shù)組中的值，來達到更行紋理存儲器中的數(shù)據(jù)的目的。這個過程中，因為數(shù)據(jù)仍然在device端進行傳遞，耗時不長。

其次，存放小波變換處理后的數(shù)據(jù)，可采用全局存儲器。全局存儲器使用的是普通顯存，整個grid的所有線程都能訪問，由于沒有片內(nèi)緩存，訪問速度慢，但通過合并訪問，可以隱藏延時。

此外，卷積過程中用到的“分解尺度函數(shù)”和“分解母函數(shù)”初始化后不再發(fā)生改變，故可將它們存放在常數(shù)存儲器中。常數(shù)存儲器大小為64 KB，是只讀的地址空間，其數(shù)據(jù)位于顯存，但擁有緩存加速。在CUDA程序中常數(shù)存儲器一般用于存儲需要頻繁訪問的只讀參數(shù)[7]。

3.3.3? 并行算法實現(xiàn)

host端控制小波變換的層數(shù)，device端執(zhí)行一次小波變換。host端關(guān)鍵偽代碼為：

for（int n=0;n

width=lWidth>>n? ? ;? //通過移位操作，控制小波變化的范圍

height=lHeight>>n ;? //通過移位操作，控制小波變化的范圍

//根據(jù)小波變換的層次，不斷調(diào)整block的數(shù)量

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

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

cudaMemcpyToArray（cu_array，0，0，d_Result，size，cudaMemcpyDeviceToDevice）; //更新紋理數(shù)據(jù)

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

cudaMemcpyToArray（cu_array，0，0，d_Result，size，cudaMemcpyDeviceToDevice）; //更新紋理數(shù)據(jù)

}

因為每次小波變換處理的區(qū)域都在發(fā)生改變，為了在device端更方便地控制線程操作，故每次循環(huán)都重新定義grid，修改grid中block的數(shù)量。

接下來，調(diào)用device端的行（列）卷積運算，進行小波變化，每完成一次行（列）卷積，都必須更新紋理數(shù)據(jù)，以供下一次循環(huán)數(shù)據(jù)調(diào)用。device端小波變換偽代碼為：

//行方向上卷積

if（x

for（int j=0; j< 某參數(shù); j++）{

m=（某計算）% width;? // width是小波變換作用的區(qū)域

tmp += tex2D（texData， （（float）m + 0.5f） ， iy） *LF[j]; //時域與分解尺度函數(shù)相乘

}

f[y*lwidth+x] = tmp; // lwidth是原圖像的寬度，f存放結(jié)果

}else if（x

for（int j=0; j< 某參數(shù); j++）{

m=（某計算）% Width; // width是小波變換作用的區(qū)域

tmp += tex2D（texData， （（float）m + 0.5f） ， iy） *HF[j]; //時域與分解母函數(shù)相乘

}

f[y*lwidth+x] = tmp; // lwidth是原圖像的寬度，f存放結(jié)果

}

上述代碼的執(zhí)行過程是并行的，即每個線程會并發(fā)地去執(zhí)行相應(yīng)的計算，然后將計算的結(jié)果保存在全局存儲器中。列方向卷積同理，不再列舉。

由上述分析可見，算法的時間復(fù)雜度由原來的四重循環(huán)變?yōu)閮芍匮h(huán)，而兩重循環(huán)的基數(shù)非常小，所以大大地縮短了算法運行時間。

4? 實驗結(jié)果分析

本程序目的在于對比小波變換在CPU和GPU上的運行效率，通過調(diào)用同一幅圖片進行小波變換，取得各自的執(zhí)行時間，最后算出加速比。

硬件上，CPU采用IntelCore i5-7200，主頻為2.5 GHz;GPU采用NVIDIA公司的GeForce GTX660，支持CUDA 4.0以上驅(qū)動版本，顯存容量2 GB。

實驗程序采用Daubechies D4小波，變換層次為3。所有數(shù)據(jù)為5次計算時間平均值，CUDA計算時間包括了數(shù)據(jù)在host端（主機）與device端（顯卡存儲器）中交換的時間。

經(jīng)過測試，小波變換在CPU和GPU上執(zhí)行的效率對比如表1所示。從表中數(shù)據(jù)可以看出，應(yīng)用CUDA技術(shù)加速小波變換，均達到兩個數(shù)量級的加速效果。

由圖3可知，隨著圖像大小的增加，CPU上小波變換所需的時間急劇上升，而GPU上執(zhí)行所需時間則增長緩慢。由圖4可見，圖像大小與加速比總體呈正向關(guān)系。

5? 結(jié)? 論

本文重點分析了小波Mallat算法并行化的可行性，并成功利用CUDA技術(shù)實現(xiàn)了高速并行小波變換算法。實驗表明，相較于CPU上的串行小波變換算法，本文算法的執(zhí)行效率提升了數(shù)十倍，最高達到了約51.6倍，且算法效率與計算量成正向關(guān)系。通過合理設(shè)計CUDA程序，將耗時的串行計算過程并行化，充分利用GPU強大的計算資源，能最大限度地縮短執(zhí)行時間，使得小波技術(shù)更好地應(yīng)用于實際工程系統(tǒng)中。

參考文獻：

[1] 陳慶奎，王海峰，那麗春，等.圖形處理器通用計算的研究綜述 [J].黑龍江大學自然科學學報，2012，29（5）：672-679.

[2] 陳茜.基于GPU的數(shù)字圖像處理并行算法的研究 [D].西安：中國科學院研究生院（西安光學精密機械研究所），2014.

[3] 白兆峰.基于GPU的JPEG2000圖像壓縮編碼技術(shù)的研究 [D].西安：中國科學院大學（中國科學院西安光學精密機械研究所），2016.

[4] 張舒，褚艷利.GPU高性能運算之CUDA [M].北京：中國水利水電出版社，2009：20-22.

[5] 陳天華.數(shù)字圖像處理 [M].北京：清華大學出版社，2007：132-133.

[6] 韓志偉，劉志剛，魯曉帆，等.基于CUDA的高速并行小波算法及其在電力系統(tǒng)諧波分析中的應(yīng)用 [J].電力自動化設(shè)備，2010，30（1）：98-101+105.

[7] 鄒巖，楊志義，張凱龍.CUDA并行程序的內(nèi)存訪問優(yōu)化技術(shù)研究 [J].計算機測量與控制，2009，17（12）：2504-2506.

作者簡介：張金霜（1987—），男，漢族，廣東茂名人，教師，畢業(yè)于華南農(nóng)業(yè)大學，碩士研究生，研究方向：計算機教育。