基于多項式運算的無累積誤差Loeffler DCT算法

沈永珞，章 媛，楊迪威，李 璇

（1.廣東財經(jīng)大學(xué) 信息學(xué)院，廣東 廣州510320；2.國家知識產(chǎn)權(quán)局專利局專利審查協(xié)作廣東中心，廣東 廣州510530；3.中國地質(zhì)大學(xué)（武漢） 數(shù)學(xué)與物理學(xué)院，湖北 武漢430074）

0 引 言

離散余弦變換 （DCT）被廣泛應(yīng)用于圖像和視頻編碼、特征提取和恢復(fù)等領(lǐng)域，并已有大量的研究通過蝶形算法及整數(shù)運算等方法，致力于提高其運算速度并降低其運算復(fù)雜度［1－4］。但是，在大數(shù)據(jù)量以及高圖像質(zhì)量要求的應(yīng)用中，比如HEVC 等，整數(shù)DCT 運算的精度影響著編解碼后視頻圖像的質(zhì)量，其原因在于運算過程中無理數(shù)參與的運算產(chǎn)生的累積誤差。因此，設(shè)計高效算法，實現(xiàn)無乘法、無運算誤差累積的DCT 運算具有極其重要的意義。

1 Loeffler DCT算法簡介

Loeffler DCT 算法最早提出于1989年，和其他快速算法相比，該算法實現(xiàn)8點一維DCT 變換只需要11次乘法和29次加法運算，其乘法器的使用數(shù)量降到了理論的極值。完整的Loeffler DCT 算法由4個階段組成，其信號運算流程如圖1所示，符號解釋見表1。其中，涉及到乘法操作的分別位于Stage－2和Stage－3中的kcn 模塊以及stage－4中的倍乘因子運算。

圖1 Loeffler DCT 算法流程

表1 符號及計算公式

對于Loeffler DCT 算法的研究一直是國內(nèi)外研究的熱點。文獻 ［5］發(fā)表了一種基于Loeffler算法的16點快速整數(shù)DCT 算法，并能取得較好的視頻編碼壓縮率。文獻［6］發(fā)表了一種改進的擴展8點Loeffler算法，通過結(jié)合量化過程中的系數(shù)，將kcn 模塊改進到蝶形運算的尾部，節(jié)省了運算量。文獻 ［7］采用CSD 編碼的方式避免Loeffler算法中的乘法操作，并設(shè)計出一種流水線結(jié)構(gòu)的IP 軟核。文獻 ［8］通過分析Blackfin533 的并行工作方式，經(jīng)過Blackfin指令集優(yōu)化，在DSP 上實現(xiàn)了Loeffler快速DCT算法的JPEG 壓縮。由于人們對圖像質(zhì)量的要求越來越高，近年來，對于精確DCT 運算的研究逐漸成為新的研究領(lǐng)域。文獻 ［9，10］發(fā)表了一種基于Arai DCT 算法的精確運算，并在硬件上設(shè)計實現(xiàn)。

已發(fā)表的對于Loeffler DCT 的研究往往集中于如何避免乘法運算，而忽略了算法中涉及的誤差，尤其是由于無理數(shù)的操作而引起的計算累積誤差。因此，本文結(jié)合現(xiàn)實應(yīng)用中對高精度和無乘法的需求，提出了一種無乘法無計算累積誤差的Loeffler DCT 算法。

2 基于多項式運算的Loeffler DCT算法實現(xiàn)

2.1 多項式運算簡介

在各種快速DCT 變換過程中，誤差來源于蝶形變換中無理數(shù)常數(shù)所參與的計算。例如，在Loeffler DCT 算法中，Stage－2階段的kcn 模塊中包含的無理數(shù)常數(shù)Cos（nπ／16）參加各種運算后，其誤差將帶入Stage－3以及Stage－4階段的計算，從而產(chǎn)生更大的累積誤差。因此，本文提出了在蝶形變換過程中，使用多項式代替上述無理數(shù)常數(shù)的方法，從而實現(xiàn)在蝶形變換計算過程中無累積誤差。并且，通過對多項式變量系數(shù)的合理選擇，將原有的復(fù)雜乘法操作化簡成為簡單的移位和加法運算。

式中：a00，a01，a10，a11為多項式的整數(shù)變量，本文中簡稱“多項式變量”，其具體數(shù)值將隨著不同的Ck而變化。對于上述公式中出現(xiàn)的Z0，Z1，Z0Z1以及后文中出現(xiàn)的Z21和Z0Z21，本文中稱為 “多項式常量”。

雖然在一定誤差范圍內(nèi)，多項式變量a00，a01，a10，a11有多種可能滿足Ck≈Pk（a00，a01，a10，a11），但是，為了計算過程中避免乘法操作，本算法中將各多項式變量的數(shù)值選擇為絕對值接近2的n 次冪的數(shù)值。因此，對于Loeffler DCT 算法中出現(xiàn)的各無理數(shù)常數(shù)，本算法中采用的其對應(yīng)的各多項式系數(shù)數(shù)值見表2。

表2 多項式變量及絕對誤差

為表達方便，算法中涉及的各無理數(shù)可用一組多項式11來表達，比如Cos （π／16）將變量 數(shù) 組 a00，a01，a10，a（））0 。通過上述將無理數(shù)轉(zhuǎn)換成多項式的表示成 0，2，－2，（表達方式后，凡是涉及到無理數(shù)Ck的運算都將轉(zhuǎn)換為多項式整數(shù)變量之間的運算。根據(jù)整數(shù)Loeffler DCT 的蝶形運算過程，本算法包含2種類型的多項式運算：

（1）A×Ck

在此類型的多項式運算中，A 為任意整數(shù)，比如經(jīng)過整數(shù)Loeffler DCT 的Stage－1階段的加減運算后的數(shù)值。由

11表示，因此，該于無理數(shù)Ck由多項式Pka00，a01，a10，a（

）類運算轉(zhuǎn)換為如下形式

其 運 算 結(jié) 果 表 示 為 （Aa00，Aa01，Aa10，Aa11）。由 于a00，a01，a10，a11的數(shù)值接近于2的n次冪，因此，可以通過簡單的移位和加減運算來代替A 與多項式整數(shù)變量之間的乘法操作，從而實現(xiàn)DCT 運算過程中的無乘法操作。

（2）M×Ci＋N×Cj

在此類型的多項式運算中，M×Ci和N×Cj的結(jié)果可經(jīng)由類型 （1）的運算后得到，分別用多項式變量數(shù)組（ai00，ai01，ai10，ai11）和 （aj00，aj01，aj10，aj11）表 示。由于上述各多項式變量對應(yīng)相同的多項式常數(shù)，因此，該類型的運算過程轉(zhuǎn)換為簡單的對應(yīng)變量部分的相加運算，其結(jié)果表示為 （ai00＋aj00，ai01＋aj01，ai10＋aj10，ai11＋aj11）。計算過程如下列公式所示

由此可見，在基于多項式的計算過程中，存在著計算累積誤差的無理數(shù)數(shù)值運算被完全避免，取而代之的是無計算累積誤差的對多項式整數(shù)變量的運算操作。因此，該算法能完全避免計算過程中的累積誤差。

2.2 改進的Loeffler DCT算法

在本文中改進的Loeffler DCT 算法中，根據(jù)運算模式的不同，可分為3種不同類型的模塊：常規(guī)數(shù)值運算模塊，多項式運算模塊和結(jié)果再生運算模塊。各計算模塊所涉及的范圍如圖2所示。

圖2 改進算法中的3個計算模塊

2.2.1 常規(guī)數(shù)值運算模塊

根據(jù)8點Loeffler DCT 算法，可以直接用常規(guī)數(shù)值運算得到y(tǒng)0和y4的結(jié)果。同時，對于蝶形運算過程中的中間結(jié)果，比如Stage－1階段運算的下半部分和Stage－2階段運算的中間部分，也可以直接通過常規(guī)數(shù)值運算得到。

2.2.2 多項式運算模塊

多項式運算模塊包含了蝶形運算中大部分的復(fù)雜運算，如Stage－2，Stage－3和Stage－4中的大部分運算。多項式運算規(guī)則如前文介紹所示，經(jīng)過多項式運算后，y1、y2、y6、y7以及和倍乘因子 槡2相乘之前的y3和y5均由多項式變量數(shù)組表達。

2.2.3 結(jié)果再生運算模塊

結(jié)果再生模塊的主要目的是將多項式運算模塊得到的結(jié)果 （a′00，a′01，a′10，a′11），計算成DCT 變換后各種后繼操作所需的具體常規(guī)數(shù)值。其計算過程為根據(jù)多項式變量數(shù)組，累加各多項式變量與其對應(yīng)多項式常數(shù)的乘積。對于y1、y2、y6和y7，其結(jié)果多項式變量數(shù)組元素分別對應(yīng)多項式常數(shù)1，Z0，Z1和Z0Z1，因此，其結(jié)果再生計算過程如式 （1）所示。

由于Loeffler DCT 算法中的倍乘因子槡2正好是本算法所選擇的多項式常數(shù)Z1，所以對于y3和y5，只需在多項式運算后得到的多項式變量數(shù)組的基礎(chǔ)上，在結(jié)果再生步驟中調(diào)整其對應(yīng)的多項式常數(shù)為Z1，Z21，Z0Z1和Z0Z21，即可實現(xiàn)與倍乘因子槡2相乘的結(jié)果。因此，y3和y5的結(jié)果再生計算過程如式 （4）所示

在結(jié)果再生運算過程中，由于大部分多項式常數(shù)仍為無理數(shù)，因此本算法考慮在一定誤差范圍內(nèi)，將各多項式常數(shù)用以2為底位權(quán)展開的表示方式參與結(jié)果再生計算，從而避免多項式變量與多項式常數(shù)的乘法操作，僅使用移位和加法。各多項式常數(shù)按位權(quán)表達及誤差見表3。

表3 多項式常數(shù)表達方式及誤差

3 運算量與實驗結(jié)果分析

3.1 運算量分析

本算法具有運算過程中無乘法操作的優(yōu)點，因此在實現(xiàn)基于此算法的硬件DCT 變換模塊能達到有較高的工作頻率，滿足圖像視頻處理中實時處理的要求。雖然多項式運算將乘法操作分解成多個移位和加法運算，但是由于各多項式變量運算的獨立性，可以通過硬件的并行處理來提高運算速度。而在硬件實現(xiàn)上，僅需要增加硬件開銷有限的整數(shù)加法器。

表4總結(jié)了本文所提出的一維8點無乘法Loeffler DCT算法的各模塊的運算開銷。從該表中可觀察到結(jié)果再生模塊所需的計算量較大，但是對于實際應(yīng)用中需要的二維8點DCT 變換，結(jié)果再生操作可以盡量避免。為減少運算開銷，行列變換的元素可以為多項式的變量數(shù)組，從而僅需使用一次結(jié)果再生操作。使用多項式運算的二維8點DCT過程如圖3所示。

表4 各模塊的運算開銷

圖3 使用多項式運算的二維8點DCT 過程

3.2 圖像質(zhì)量分析

為了驗證本文提出的算法對圖像編解碼后圖像質(zhì)量的影響，在實驗過程中，首先根據(jù)本文提出的算法將圖像進行整數(shù)DCT 處理，然后使用浮點運算的Loeffler IDCT 進行圖像恢復(fù)，最后將重建圖像與原圖像進行比較，以峰值信噪比PSNR 作為質(zhì)量的量化衡量標準。其中，整數(shù)運算過程中，使用長度為12 位的整數(shù)數(shù)據(jù) （包含2位小數(shù)部分）。實驗采用的圖像如圖4 所示 （圖像大小256×256），其結(jié)果PSNR值見表5。實驗結(jié)果表明，該算法能在避免使用乘法操作的同時，圖像質(zhì)量也能得到較高的保證。

圖4 測試圖片

4 結(jié)束語

本文為了解決原始Loeffler DCT 算法中存在乘法操作和對無理數(shù)操作產(chǎn)生累積誤差的問題，提出了一種基于多項式計算的DCT 算法。通過使用多項式參數(shù)數(shù)組表示無理數(shù)的方法，將涉及到無理數(shù)的運算轉(zhuǎn)換成簡單的多項式變量之間的運算，從而避免了蝶形運算過程中計算誤差的累積。同時，通過選擇合適的多項式變量數(shù)值，在多項式運算過程中避免了復(fù)雜的乘法運算。實驗表明，在增加有限的加法和移位操作的前提下，該算法能保證較好的圖像質(zhì)量。

［1］Lee DU，Kim H，Rahimi M，et al.Enginery－efficient image compression for resource－constrained platforms ［J］.IEEE Trans Image Process，2009，18 （9）：2100－2113.

［2］YAO Qingdong，BI Houjie，WANG Zhaohua，et al.Image coding primer ［M］.3rd ed.Beijing：Tsinghua University Press，2006 （in Chinese）. ［姚慶東，畢厚杰，王兆華，等.圖像編碼基礎(chǔ) ［M］.3版.北京：清華大學(xué)出版社，2006.］

［3］LI Kangshun，WEI Yunshan，ZHANG Wensheng.Wraped DCT image compression based on evolutionary algorithm ［J］.Journal of Image and Graphics，2013，18 （2）：169－175 （in Chinese）. ［李康順，韋蘊珊，張文生.基于演化算法的卷曲DCT 圖像壓縮［J］.中國圖象圖形學(xué)報，2013，18 （2）：169－175.］

［4］PANG Rong，LIU Yu，HOU Zhengxin，WANG Shaochu.Image coding and reconstruction via compressed sensing based on partial DCT coefficients ［J］.ACTA Automatic Sinica，2011，37 （6）：674－681 （in Chinese）. ［潘榕，劉昱，侯正信，汪少初.基于局部DCT 系數(shù)的圖像壓縮感知編碼與重構(gòu)［J］.自動化學(xué)報，2011，37 （6）：674－681.］

［5］Fong CK，Cham WK.LLM integer cosine transform and its fast algorithm ［J］.IEEE Trans on Circus and System for Video Technology，2012，22 （6）：844－854.

［6］Wu ZG，Sha J，Wang ZF，et al.An improved scaled DCT architecture［J］.IEEE Trans on Consumer Electronics，2009，55 （2）：685－689.

［7］GUO Baozeng，NIU Li，LIU Zhiming.Dsign of 2－D DCT IP soft core based on Loeffler’s algorithm ［J］.Mircroelectronics＆Computer，2011，28 （2）：136－144 （in Chinese）. ［郭寶增，牛力，劉志明.基于Loeffler算法的2－D DCT IP 軟核設(shè)計 ［J］.微電子學(xué)與計算機，2011，28 （2）：136－144.］

［8］LIU Xuan，LIU Ziyi，XU Kaoji.JPEG compression using Loeffler 1D－DCT algorithm based on DSP ［J］.Journal of Data Acquisition ＆ Processing，2009，24 （S1）：61－64 （in Chinese）.［劉軒，劉子軼，徐考基.基于Loeffler快速DCT 算法的JPEG 壓縮在DSP 上的實現(xiàn) ［J］.數(shù)據(jù)采集與處理，2009，24 （S1）：61－64.］

［9］Edirisuriya A，Madanayake A，Dimitrove V，et al.VLSI architecture for 8－point AI－based Arai DCT having：Low areatime complexity and power at improved accuracy ［J］.Journal of Low Power Electronics and Applications，2012，2 （2）：127－142.

［10］Madanayake A，Cintra RJ，Onen D，et al.A Row－parallel 88 2－D DCT architecture using algebraic integer－based exact computation ［J］.IEEE Trans on Circuits and System for Video Technology，2012，22 （6）：915－929.