基于CUDA技術(shù)的DCT并行算法研究與實(shí)現(xiàn)

張金霜 黃旭彬

摘要：JPEG有損壓縮算法以DCT變換算法為核心，但DCT變換算法較為耗時(shí)，為提升圖像壓縮效率，提出利用基于GPU平臺(tái)的CUDA技術(shù)對(duì)DCT算法做并行優(yōu)化。通過(guò)分析DCT變換算法的原理，發(fā)現(xiàn)DCT算法具有很好的并行性，優(yōu)化空間較大，于是利用CUDA技術(shù)實(shí)現(xiàn)高速DCT并行算法。實(shí)驗(yàn)結(jié)果表明，在一張2592x2592分辨率的圖像做DCT變換，加速比能達(dá)到124.5，算法效率得到明顯提升，且壓縮效果無(wú)顯著差異。

關(guān)鍵詞：CUDA;GPU;并行程序設(shè)計(jì);DCT算法;圖像壓縮

中圖分類號(hào)：TP393 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2020）33-0198-04

開(kāi)放科學(xué)（資源服務(wù)）標(biāo)識(shí)碼（OSID）：

1 背景

離散余弦變換（ Discrete Cosine Transform，DCT）類似于傅立葉變換，由Ahmed在1974年首次提出。離散余弦變換經(jīng)常被用作信號(hào)和圖像處理，如對(duì)信號(hào)和圖像進(jìn)行有損數(shù)據(jù)壓縮。在JPEG圖像壓縮算法中，有一種是以離散余弦變換為基礎(chǔ)的有損壓縮算法[1]。作為JPEG壓縮算法的核心，DCT是其中最耗時(shí)的部分，且占有較大的工作量比重。

目前GPU作為一種擅長(zhǎng)處理大規(guī)模并行數(shù)據(jù)計(jì)算的設(shè)備，已受到越來(lái)越多人的關(guān)注。相較于CPU，GPU有著強(qiáng)大的計(jì)算能力和超高度存儲(chǔ)器帶寬，可以使用數(shù)量較多的執(zhí)行單元，實(shí)現(xiàn)數(shù)據(jù)間的并發(fā)計(jì)算。而CUDA技術(shù)能充分利用GPU的并行計(jì)算能力，且降低了開(kāi)發(fā)人員的使用門檻，基于GPU的CUDA并行加速技術(shù)也廣泛應(yīng)用于圖像處理領(lǐng)域[2-4]。因此，本文將DCT串行算法改為基于CPU/GPU的異構(gòu)算法，通過(guò)GPU內(nèi)部大量獨(dú)立的計(jì)算單元，實(shí)現(xiàn)計(jì)算加速。

2 CUDA體系結(jié)構(gòu)

CUDA是由顯卡廠商N(yùn)VIDIA推出的通用并行計(jì)算框架，該架構(gòu)能使有效調(diào)用GPU資源解決復(fù)雜的計(jì)算問(wèn)題。CUDA允許開(kāi)發(fā)人員使用C/C++語(yǔ)言和CUDA擴(kuò)展庫(kù)的形式編寫(xiě)程序，避免了傳統(tǒng)方法中將通用計(jì)算轉(zhuǎn)換為圖形處理的過(guò)程，大大降低程序設(shè)計(jì)的難度和系統(tǒng)維護(hù)成本。

CUDA實(shí)現(xiàn)了CPU和GPU的混合架構(gòu)編程。在CUDA架構(gòu)下，程序被分成兩個(gè)部分：host端程序和device端程序。host端程序運(yùn)行在CPU上;device端程序又叫“kernel”，運(yùn)行在GPU上。通常先由host端程序準(zhǔn)備數(shù)據(jù)，對(duì)設(shè)備進(jìn)行必要的初始化，把數(shù)據(jù)從內(nèi)存?zhèn)魉偷斤@存中;執(zhí)行device端程序，待GPU運(yùn)行完后，host端程序再?gòu)娘@存中把結(jié)果取回內(nèi)存中，過(guò)程如圖1所示[5]。

CUDA的線程模型分三個(gè)層次：thread、block、grid。其中，kernel以線程網(wǎng)格（grid）的形式組織，每個(gè)線程網(wǎng)格由若干個(gè)線程塊（block）組成，而每個(gè)線程塊又由若干個(gè)線程（thread）組成。

block內(nèi)的thread之間是并行執(zhí)行的，grid內(nèi)的block之間也是并行執(zhí)行的，而grid與grid之間是串行執(zhí)行的，即GPU同時(shí)只能處理一個(gè)grid。grid包含一組block，grid是kernel在GPU上的一次執(zhí)行，也就是說(shuō)，一個(gè)kernel對(duì)應(yīng)著一個(gè)grid。grid內(nèi)部的block可以按照一維或者二維來(lái)組織，block之間共享全局內(nèi)存。

3 DCT算法理論

二維頻率與矩陣c（u，v）的左上角代表二維正交展開(kāi)中的低階項(xiàng)，它代表圖像的低頻分量;c（u，v）的右下角代表二維正交展開(kāi)中的高階項(xiàng)，它代表圖像的高頻分量。圖像的低頻分量，反映圖像灰度的緩慢變化;相反，圖像的高頻分量反映圖像灰度的跳變，又稱輪廓或邊緣[6]。

經(jīng)過(guò)DCT變換后的矩陣數(shù)據(jù)自然數(shù)為頻率系數(shù)，這些系數(shù)中，位于左上角的F（0，0）的值最大，稱為直流系數(shù)（DC coeffi-cient）;其余的63個(gè)頻率系數(shù)多為一些接近于0的正負(fù)浮點(diǎn)數(shù)，稱為交流系數(shù)（AC coeffcients）[7]。利用這一特點(diǎn)就可以實(shí)現(xiàn)圖像數(shù)據(jù)壓縮，對(duì)其進(jìn)行量化，僅僅保留低頻分量的左上角，其余分量拋棄，這樣就可以達(dá)到圖像數(shù)據(jù)壓縮的目的。在反變換的時(shí)候，拋棄的點(diǎn)用零填充。

4 基于CUDA實(shí)現(xiàn)DCT算法加速

4.1可行性分析

CUDA架構(gòu)的優(yōu)勢(shì)在于處理數(shù)據(jù)并行的問(wèn)題。通過(guò)生成盡可能多的線程來(lái)處理，每一個(gè)線程的處理是基本相同的，而處理的數(shù)據(jù)是可以并行的。符合這種特點(diǎn)的算法，能在CUDA架構(gòu)下最大化并行執(zhí)行。

根據(jù)JPEG標(biāo)準(zhǔn)，圖像被分成一個(gè)個(gè)8x8的像素塊進(jìn)行DCT變換。針對(duì)每一個(gè)8x8的像素塊，它們的處理是相同的，都是進(jìn)行DCT變換。而每一個(gè)8x8的像素塊是相互獨(dú)立的，塊與塊之間沒(méi)有關(guān)聯(lián)。而變換矩陣是確定的，與變換對(duì)象無(wú)關(guān)。因此對(duì)各8x8像素塊進(jìn)行DCT變換是數(shù)據(jù)并行的。

4.2 算法總體流程

1）算法從host端開(kāi)始，首先在內(nèi)存中分配兩個(gè)具有與被處理圖像像素?cái)?shù)據(jù)相同大小的空間，分別用于存放被處理圖像的源數(shù)據(jù)和處理結(jié)果，并把圖像像素?cái)?shù)據(jù)讀人到內(nèi)存中。圖像要在GPU中進(jìn)行處理，于是在顯存中申請(qǐng)兩個(gè)具有與被處理圖像像素?cái)?shù)據(jù)相同大小的空間，分別用于存放圖像源數(shù)據(jù)和處理結(jié)果，并一次性把圖像像素?cái)?shù)據(jù)從內(nèi)存復(fù)制到顯存中。

2）把每一個(gè)8x8塊讀人到共享存儲(chǔ)器中。采用多線程讀人，每個(gè)線程讀人一個(gè)像素。讀人完成后，進(jìn)行DCT變換的計(jì)算。采用多線程進(jìn)行計(jì)算，每個(gè)線程計(jì)算一個(gè)系數(shù)。

3）DCT變換完成后，把結(jié)果寫(xiě)入到全局存儲(chǔ)器中，同樣是每個(gè)線程寫(xiě)入一個(gè)系數(shù)。至此device端的工作完成。

4）host端再把結(jié)果數(shù)據(jù)從顯存復(fù)制到內(nèi)存中，DCT變換算法完成。

4.3 算法詳細(xì)設(shè)計(jì)

1）線程劃分

把整個(gè)圖像作為一個(gè)grid，由于對(duì)圖像采取的是8x8塊的DCT變換，塊與塊之間可以獨(dú)立進(jìn)行處理，不同塊沒(méi)有關(guān)聯(lián)。因此，這里采取每個(gè)8x8塊用一個(gè)block來(lái)處理，而block中包含8x8個(gè)thread（共64個(gè)thread），每個(gè)thread處理一個(gè)像素的計(jì)算。工作過(guò)程如圖2所示。

2）存儲(chǔ)器分配

根據(jù)各存儲(chǔ)器的特點(diǎn)，先把圖像數(shù)據(jù)載人到全局存儲(chǔ)器中。由于全局存儲(chǔ)器的訪問(wèn)速度慢，而且沒(méi)有緩沖，每一個(gè)block內(nèi)，各thread要取的數(shù)據(jù)有重復(fù)之處，如果線程直接從全局存儲(chǔ)器取數(shù)據(jù)，將會(huì)遭遇性能瓶頸。因此，由各thread從全局存儲(chǔ)器中讀人到共享存儲(chǔ)器中，各自線程讀取自身坐標(biāo)對(duì)應(yīng)的像素值。在各thread讀入像素值后，要進(jìn)行柵障同步，以保證同一個(gè)block內(nèi)所有像素都已經(jīng)讀入，這樣才能進(jìn)行下一步的DCT變換。

為了進(jìn)一步優(yōu)化，還可以把全局存儲(chǔ)器的中數(shù)據(jù)映射為紋理，紋理有片內(nèi)緩沖，而且能進(jìn)行二維尋址，能簡(jiǎn)化DCT算法的尋址，提高訪存的速度。同時(shí)，紋理數(shù)據(jù)按16的倍數(shù)對(duì)齊，以符合“合并訪問(wèn)”的要求，進(jìn)一步提高訪存效率。

此外，DCT變換所需的變換矩陣A要先載人到常數(shù)存儲(chǔ)器，常數(shù)存儲(chǔ)器同樣具有片內(nèi)緩沖。在代碼中，用一constant一聲明限定符直接聲明為常數(shù)存儲(chǔ)器變量即可[8]。

3）并行DCT算法實(shí)現(xiàn)

根據(jù)等式（5），DCT變換要計(jì)算兩次矩陣乘法，分別是ATXArX和（ArX）A。矩陣相乘的算法，可以通過(guò)線程的索引來(lái)確定其計(jì)算的數(shù)據(jù)地址。在CUDA中，線程是能獲得自身在block中的索引的，基于這一特點(diǎn)，在矩陣乘法運(yùn)算中，可以很方便地在共享存儲(chǔ)器中尋址。每個(gè)矩陣都需要用到轉(zhuǎn)置矩陣A，而矩陣A被載入到常數(shù)存儲(chǔ)器，GPU能對(duì)常數(shù)存儲(chǔ)器進(jìn)行片內(nèi)緩沖，訪問(wèn)效率很高。

以下為DCT變換中兩次求矩陣乘積的核心代碼：

//計(jì)算ATX，結(jié)果保存到CurBlockLoca12

float curelem=0;//臨時(shí)保存當(dāng)前元素的值

int DCTv8matrixlndex=ty;//轉(zhuǎn)置矩陣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;/，保存結(jié)果

_syncthreads0;//等待矩陣的所有元素計(jì)算完畢

//計(jì)算（ATX） A（ATX）A，結(jié)果保存到CurBlockLocall

curelem=0;//臨時(shí)保存當(dāng)前元素的值

int CurBlockLoca121ndex

=

（ty

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

DCTv8matrixlndex=tx;//A矩陣索引

for （int i=0; i

（

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

CurBlockLocal21ndex+=1;//取行，因此步長(zhǎng)為1

DCTv8matrixlndex+=BLOCK_SIZE;，/取列，因此步長(zhǎng)為BLOCK__ SIZE

}

CurBlockLocall[（ty<

_syncthreads0;//柵障同步，確保所有系數(shù)計(jì)算完成

5 實(shí)驗(yàn)結(jié)果分析

5.1 數(shù)據(jù)分析

為了測(cè)試應(yīng)用CUDA技術(shù)后的運(yùn)行效率，本文編寫(xiě)了兩個(gè)版本的程序，分別在CPU運(yùn)行（單線程），和在GPU上運(yùn)行（基于CUDA技術(shù)）。通過(guò)分別調(diào)用兩個(gè)程序，針對(duì)同一幅圖片進(jìn)行DCT變換，并計(jì)算其執(zhí)行時(shí)間。

在計(jì)算執(zhí)行時(shí)間上，本文采用CUDA的cutStartTimer（tim-er）、cutStopTimer（timer）和cutGetTimerValue（timer）函數(shù)，分別表示開(kāi)始計(jì)時(shí)、停止計(jì)時(shí)和讀取時(shí)間，計(jì)時(shí)精度小于Ims。

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

經(jīng)過(guò)測(cè)試，DCT算法在CPU和GPU上執(zhí)行的效率對(duì)比如表1所示。從表中數(shù)據(jù)可以看出，應(yīng)用CUDA技術(shù)加速DCT變換，加速比最高能達(dá)到124.5。

5.2 程序演示

本文通過(guò)模塊化的方式，把DCT變換算法這個(gè)功能作為一個(gè)獨(dú)立模塊。DCT變換算法的CUDA實(shí)現(xiàn)版本和CPU實(shí)現(xiàn)版本（C語(yǔ)言洛編譯成一個(gè)DLL文件。主程序模塊則由C#語(yǔ)言編寫(xiě)，編譯成一個(gè)可執(zhí)行文件。

以下程序分別演示了應(yīng)用CUDA技術(shù)在GPU上進(jìn)行圖像壓縮，和同樣的算法在CPU上進(jìn)行計(jì)算的效果。左邊的圖是原圖，中間的圖是經(jīng)過(guò)DCT變換后的系數(shù)，右邊的圖是經(jīng)過(guò)DCT——量化——逆量化-IDCT后的結(jié)果。

6 結(jié)束語(yǔ)

本文重點(diǎn)分析了DCT變換算法的原理，得出DCT變換算法的矩陣式表示，并成功利用CUDA技術(shù)實(shí)現(xiàn)了高速DCT變換算法。實(shí)驗(yàn)表明，相較于CPU上的串行DCT變換算法，本文算法的執(zhí)行效率得到了很大提升。

本文也提供了一種程序優(yōu)化思路——對(duì)于耗時(shí)的串行程序，只要滿足并行性的條件，都可嘗試采用CUDA技術(shù)進(jìn)行改進(jìn)。而把CUDA部分的程序編譯成DLL文件，可靈活地被第三方項(xiàng)目調(diào)用，使得CUDA技術(shù)適用于更多技術(shù)場(chǎng)景。

參考文獻(xiàn)：

[1]萬(wàn)豐，苑豪杰，宮威，基于離散余弦變換的數(shù)字圖像壓縮技術(shù)研究[J].自動(dòng)化應(yīng)用，2020（3）：65-67.

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

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

[4]曾浩洋.圖像處理的GPU加速技術(shù)的研究及實(shí)現(xiàn)[D].綿陽(yáng)：西南科技大學(xué)，2019.

[5]張舒，褚艷利.GPU高性能運(yùn)算之CUDA[M].北京：中國(guó)水利水電出版社，2009.

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

[7]佘秋菊.基于DCT變換的JPEG圖像壓縮及其MATLAB實(shí)現(xiàn)[J].科技信息（學(xué)術(shù)研究），2008（36）：557-558.

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

【通聯(lián)編輯：謝媛媛】

作者簡(jiǎn)介：張金霜（1987-），男，廣東茂名人，助教，碩士，主要研究方向?yàn)榫W(wǎng)絡(luò)安全、智能優(yōu)化算法;黃旭彬（1981-），廣東高州人，講師，碩士，主要研究方向?yàn)橛?jì)算機(jī)應(yīng)用技術(shù)。