C語言實(shí)現(xiàn)DCT變換編碼1

第.4:“Wherearithmeticprecisionisnotspecified,suchasthecalculationoftheIDCT,theprecisionshallbesufficientsothatsignificanterrorsdonotoccurinthefinalintegervalues.”逆DCT變換的過程這里不再詳述，需要實(shí)現(xiàn)這個(gè)的可以去參考這個(gè)標(biāo)準(zhǔn)。在實(shí)際應(yīng)用中一般通過兩次1-DIDCT變換來完成2-DIDCT變換，這種方法通常被稱為行－列法。一般來說，后者在結(jié)構(gòu)上的對稱性更好，并且可以重復(fù)使用硬件資源，所以在我們的芯片設(shè)計(jì)選用一種行-列法來進(jìn)行IDCT單元的結(jié)構(gòu)研究。二維IDCT可以分解成二次一維IDCT運(yùn)算，如以下公式。在結(jié)構(gòu)上，上兩式所定義的運(yùn)算使用了相同的運(yùn)算“核”(如以下公式所示),它們具有相似性。因此利用三角函數(shù)的各種關(guān)系，可以得到“核”的快速算法。其中，為了便于理解，可將快速算法表示成蝶形圖，如下圖。將對模塊進(jìn)行一維IDCT變換的結(jié)果存儲(chǔ)起來，轉(zhuǎn)置輸出，再進(jìn)行一次IDCT變換，即為相應(yīng)的二維IDCT變換。圖4－8折疊結(jié)構(gòu)的二維IDCT單元一行數(shù)據(jù)(一行有8個(gè)像素?cái)?shù)據(jù))在該單元中的處理流程是：１—>２—>３—>４—>５—>６—>７—>８。DCT變換探究1前言

此文適合于那些對DCT或?qū)aar小波的Mallat算法有一定了解的人。

由于我還是高一新丁，文學(xué)底子很薄弱，對于一些技術(shù)方面的知識(shí)，我是有口說不出，無法用文字表達(dá)出來，因此這里提供的知識(shí)只是我所知道的1/4左右，還有3/4我不知該如何表達(dá)，特別是第三節(jié)“深入研究DCT”，我個(gè)人認(rèn)為簡直是淺入！

如果你只是菜鳥，不但想看懂此文，而且還要看懂其他的類似文章，那么我教你一個(gè)最快的學(xué)習(xí)方法：

設(shè)X={10,20}

分解的方法：低頻=10+20=30，高頻=10-20=-10，

即Y={30,-10}

合并的方法：X(0)=(低頻+高頻)/2=(30+(-10))/2=10,X(1)=X(0)-高頻=10-(-10)=20

即X={10,20}

只要搞清楚低頻和高頻是怎么來的和如何合并的即可。2DCT簡介

DCT全名為DiscreteCosineTransform，中文名為離散余弦變換。在眾人皆知的JPEG編碼中，就是使用了DCT來壓縮圖像的。為什么DCT可以壓縮圖像？我想這個(gè)問題有很多人都想知道，但其實(shí)這是錯(cuò)誤的說法！因?yàn)镈CT在圖像壓縮中僅僅起到扶助的作用，給它n個(gè)數(shù)據(jù)，經(jīng)變換后仍然會(huì)得出n個(gè)數(shù)據(jù)，DCT只不過消除了這n個(gè)數(shù)據(jù)的冗余性和相關(guān)性。

即，用很少的數(shù)據(jù)就能大致還原出這n個(gè)數(shù)據(jù)，其他的一些DCT系數(shù)只起到修正的作用，可有可無。

DCT有一個(gè)缺點(diǎn)，就是計(jì)算量很大！因?yàn)槿绻凑誅CT的標(biāo)準(zhǔn)變換公式(二維)來實(shí)現(xiàn)8x8點(diǎn)陣的變換需要將近上萬次計(jì)算！后來提出了一種優(yōu)化方法，即將二維的DCT分解為兩個(gè)一維的DCT，這樣一來計(jì)算量就可以減少為原來的1/4。但是計(jì)算量依然巨大，不具有使用價(jià)值，后來在1988年有人提出了一種快速算法叫AAN,它也是將二維的DCT分解成一維的形式，但是二維計(jì)算量已減少到只有600來次了，JPG和MPEG編碼中的DCT就是使用AAN算法實(shí)現(xiàn)的。

DCT還有一個(gè)缺點(diǎn)，就是不能無損變換，因?yàn)镈CT系數(shù)都是一些無理數(shù)，目前為止，依然無法解決。3深入研究首先讓我們來看看AAN算法的第一階級(jí)變換代碼：

ForI=0To3

J=7-I

Y(I)=X(I)+X(J)

Y(J)=X(I)-X(J)

NextI

設(shè)X={10,20,30,40,50,60,70,80}

那么Y={90,90,90,90,-10,-30,-50,-70}

可以看出，這一階級(jí)的低頻部分(相加得出的數(shù)據(jù))全部相等，而高頻部分則呈線性或者是有規(guī)律的。DCT之所以能以較少的數(shù)據(jù)大致還原圖像，就是因?yàn)橥ㄟ^預(yù)測高頻部分而達(dá)到的。那么為何高頻部分可以預(yù)測呢？請仔細(xì)看上面的代碼，可以看出DCT是由外到內(nèi)來進(jìn)行處理的，由于像素及像素間有一定的關(guān)聯(lián)性，所以靠的越近的像素之間的差就應(yīng)該越小，越遠(yuǎn)就因該越大，但也并不是說所有的數(shù)據(jù)都具有這種規(guī)律，因此DCT預(yù)測出來的高頻數(shù)據(jù)就會(huì)和原高頻數(shù)據(jù)不大相同，它們之間的差便是第二節(jié)提出的修正數(shù)據(jù)。第二階級(jí)變換則是在第一階級(jí)變換的基礎(chǔ)上再次分解出低、高頻，和預(yù)測高頻，得出修正值。第三階級(jí)……。最后，再將DCT系數(shù)按照重要程度由大到小，由左到右，重排列即可。

例：X={10,20,30,40,50,60,70,80}

經(jīng)過FDCT后：Y={127,-64,0,-7,-0,-2,0,-1}

其中127是最最重要的，而-64次之，以此類推?？梢园l(fā)現(xiàn),-7,-2,-1的能量都很小，說明這三個(gè)修正值可以忽略，當(dāng)忽略后，

得Y={127,-64,0,0,0,0,0,0}

經(jīng)過IDCT后：X={14,18,27,39,51,63,72,76}

這及原始數(shù)據(jù)：X={10,20,30,40,50,60,70,80}

是非常接近的，肉眼很難發(fā)覺。4為何JPEG2000放棄DCT

在JPEG2000里，放棄了基于塊的DCT，而改為了小波變換。為何要放棄DCT呢？我認(rèn)為最根本的原因還是跟DCT的計(jì)算量有關(guān)，第二節(jié)已經(jīng)指出，為了減少計(jì)算量，我們不得不使用只能處理8X8點(diǎn)陣的AAN快速算法（目前，也只有基于8X8點(diǎn)陣的），對于一幅圖像，必須將其分割成無數(shù)個(gè)8X8大小的“塊”，對塊進(jìn)行變換。在低碼率下，就會(huì)產(chǎn)生方塊效應(yīng)，要解決這個(gè)問題，唯有不使用基于區(qū)塊的AAN快速算法，而是使用直接變換法，但計(jì)算量驚人！由于小波變換計(jì)算量很少，便于直接處理圖像數(shù)據(jù)，因此就不會(huì)產(chǎn)生塊效應(yīng)，但假如用小波也進(jìn)行基于8X8點(diǎn)陣的塊變換，在低碼率下，同樣也會(huì)有塊效應(yīng)！只要是基于塊變換的，那么在低碼率下就會(huì)出現(xiàn)塊效應(yīng)，無論是DCT還是小波。因此，如果忽略DCT直接處理的計(jì)算量問題的話，我認(rèn)為壓縮效率會(huì)比JPEG2000更好?。ň唧w原因暫不討論）5DCT的改進(jìn)

下面的代碼是我對DCT變換的改進(jìn)，它具有以下特性

l無損變換

l計(jì)算量少

l原位計(jì)算

經(jīng)改進(jìn)后，它已不再叫作DCT了，可以認(rèn)為是一種新的算法，只不過是在DCT的基礎(chǔ)上修改而來。

以下是正變換：X為輸入端，Y為輸出端

設(shè)X={10,20,30,40,50,60,70,80}那么Y={45,40,0,0,0,0,0,0}

'第一階級(jí)

ForI=0To3

J=7-I

Y(J)=X(J)-X(I)

Y(I)=X(I)+Fix(Y(J)/2)

NextI

'第二階級(jí)

ForH=0To4Step4

ForI=0To1

J=3-I

X(J+H)=Y(J+H)-Y(I+H)

X(I+H)=Y(I+H)+Fix(X(J+H)/2)

NextI

NextH

'第三階級(jí)

ForI=0To6Step2

Y(I+1)=X(I+1)-X(I)

Y(I)=X(I)+Fix(Y(I+1)/2)

NextI

'預(yù)測

Y(3)=Y(3)-Y(2)

Y(6)=Y(6)-Y(7)

Y(7)=Y(7)-Y(4)

'重要性排序及AAN一樣，皆為{0,4,2,6,1,5,7,3}，此略為何能無損？為何能原位？和具體實(shí)現(xiàn)原理暫時(shí)略，以后我會(huì)補(bǔ)上6參考

[1]丁貴廣，計(jì)文平，郭寶龍VisualC++6.0數(shù)字圖像編碼p44,p57,p170快速DCT變換仿效FFT的FDCT方法有及DCT無關(guān)的復(fù)數(shù)運(yùn)算部分，選用代數(shù)分解法可以降低運(yùn)算量，達(dá)到高速運(yùn)算的目的。代數(shù)分解法實(shí)現(xiàn)如下：對一維DCT表達(dá)式直接展開，尋找各點(diǎn)表達(dá)式中共同項(xiàng)，仿FFT蝶形關(guān)系，將表達(dá)式中的共同項(xiàng)作為下一級(jí)節(jié)點(diǎn)，依次進(jìn)行多次，最后得到變換結(jié)果。一、DCT部分例子：Definecos(n*pi/16)Cn

F(0，v)=0.5*C(0)*[x(0)+x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)]F(1，v)=0.5*C(0)*[x(0)*C1+x(1)*C3+x(2)*C5+x(3)*C7+x(4)*C9

+x(5)*C11+x(6)*C13+x(7)*C15]

=0.5*{[x(0)-X(7)]C1+[X(1)-X(6)]*C3+[X(2)－x(5)]*C5

+[x(3)-x(4)]*C7]從上面的式子可以看到07,16,25,34可以作為第一次運(yùn)算的相加節(jié)點(diǎn)，將所有節(jié)點(diǎn)的表達(dá)式列出后，可發(fā)現(xiàn)一個(gè)規(guī)律，得到一蝶形圖，按之編程，如下：#define

C10.9808

#define

C20.9239

#define

C30.8315

#define

C40.7071

#define

C50.5556

#define

C60.3827

#define

C70.1951//先做行DCT

voidfdctrow(double*blk)

{

doubleS07,S16,S25,S34,S0734,S1625;

doubleD07,D16,D25,D34,D0734,D1625;S07=blk[0]+blk[7];

S16=blk[1]+blk[6];

S25=blk[2]+blk[5];

S34=blk[3]+blk[4];

S0734=S07+S34;

S1625=S16+S25;D07=blk[0]-blk[7];

D16=blk[1]-blk[6];

D25=blk[2]-blk[5];

D34=blk[3]-blk[4];

D0734=S07-S34;

D1625=S16-S25;blk[0]=0.5*(C4*(S0734+S1625));

blk[1]=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);

blk[2]=0.5*(C2*D0734+C6*D1625);

blk[3]=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);

blk[4]=0.5*(C4*(S0734-S1625));

blk[5]=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);

blk[6]=0.5*(C6*D0734-C2*D1625);

blk[7]=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);//再做列DCTvoidfdctcol(double*blk)

{

doubleS07,S16,S25,S34,S0734,S1625;

doubleD07,D16,D25,D34,D0734,D1625;S07=blk[0*8]+blk[7*8];

S16=blk[1*8]+blk[6*8];

S25=blk[2*8]+blk[5*8];

S34=blk[3*8]+blk[4*8];

S0734=S07+S34;

S1625=S16+S25;D07=blk[0*8]-blk[7*8];

D16=blk[1*8]-blk[6*8];

D25=blk[2*8]-blk[5*8];

D34=blk[3*8]-blk[4*8];

D0734=S07-S34;

D1625=S16-S25;blk[0*8]=0.5*(C4*(S0734+S1625));

blk[1*8]=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);

blk[2*8]=0.5*(C2*D0734+C6*D1625);

blk[3*8]=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);

blk[4*8]=0.5*(C4*(S0734-S1625));

blk[5*8]=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);

blk[6*8]=0.5*(C6*D0734-C2*D1625);

blk[7*8]=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);voidfdct(double*block)

{

inti;

for(i=0;i<8;i++)

fdctrow(block+8*i);

for(i=0;i<8;i++)

fdctcol(block+i);

}二、IDCT部分

圖片來源：G.G.Pechanek,C.W.Kurak,C.J.Glossner,C.H.L.Moller,andS.J.WalshIBMMicroelectronicsDivision,ResearchTrianglePark,N.C."M.F.A.S.T.:AHIGHLYPARALLELSINGLECHIPDSPWITHA2DIDCTEXAMPLE"圖中未給出系數(shù)，需自行算出，仿上述DCT方法直接展開表達(dá)式搜尋規(guī)律即可。編程如下：#define

C10.9808

#define

C20.9239

#define

C30.8315

#define

C40.7071

#define

C50.5556

#define

C60.3827

#define

C70.1951//對行做DCTvoididctrow(double*blk)

{

doubletmp[16];

//firststep

tmp[0]=blk[0]*C4+blk[2]*C2;

tmp[1]=blk[4]*C4+blk[6]*C6;

tmp[2]=blk[0]*C4+blk[2]*C6;

tmp[3]=-blk[4]*C4-blk[6]*C2;

tmp[4]=blk[0]*C4-blk[2]*C6;

tmp[5]=-blk[4]*C4+blk[6]*C2;

tmp[6]=blk[0]*C4-blk[2]*C2;

tmp[7]=blk[4]*C4-blk[6]*C6;

tmp[8]=blk[1]*C7-blk[3]*C5;

tmp[9]=blk[5]*C3-blk[7]*C1;

tmp[10]=blk[1]*C5-blk[3]*C1;

tmp[11]=blk[5]*C7+blk[7]*C3;

tmp[12]=blk[1]*C3-blk[3]*C7;

tmp[13]=-blk[5]*C1-blk[7]*C5;

tmp[14]=blk[1]*C1+blk[3]*C3;

tmp[15]=blk[5]*C5+blk[7]*C7;

//secondstep

tmp[0]=0.5*(tmp[0]+tmp[1]);

tmp[1]=0.5*(tmp[2]+tmp[3]);

tmp[2]=0.5*(tmp[4]+tmp[5]);

tmp[3]=0.5*(tmp[6]+tmp[7]);

tmp[4]=0.5*(tmp[8]+tmp[9]);

tmp[5]=0.5*(tmp[10]+tmp[11]);

tmp[6]=0.5*(tmp[12]+tmp[13]);

tmp[7]=0.5*(tmp[14]+tmp[15]);

//thirdstep

blk[0]=tmp[0]+tmp[7];

blk[1]=tmp[1]+tmp[6];

blk[2]=tmp[2]+tmp[5];

blk[3]=tmp[3]+tmp[4];

blk[4]=tmp[3]-tmp[4];

blk[5]=tmp[2]-tmp[5];

blk[6]=tmp[1]-tmp[6];

blk[7]=tmp[0]-tmp[7];

/*

blk[0]=0.5*(Y0C4+Y2C2+Y4C4+Y6C6+Y1C1+Y3C3+Y5C5+Y7C7);

blk[1]=0.5*(Y0C4+Y2C6-Y4C4-Y6C2+Y1C3-Y3C7-Y5C1-Y7C5);

blk[2]=0.5*(Y0C4-Y2C6-Y4C4+Y6C2+Y1C5-Y3C1+Y5C7+Y7C3);

blk[3]=0.5*(Y0C4-Y2C2+Y4C4-Y6C6+Y1C7-Y3C5+Y5C3-Y7C1);

blk[4]=0.5*(Y0C4-Y2C2+Y4C4-Y6C6-Y1C7+Y3C5-Y5C3+Y7C1);

blk[5]=0.5*(Y0C4-Y2C6-Y4C4+Y6C2-Y1C5+Y3C1-Y5C7-Y7C3);

blk[6]=0.5*(Y0C4+Y2C6-Y4C4-Y6C2-Y1C3+Y3C7+Y5C1+Y7C5);

blk[7]=0.5*(Y0C4+Y2C2+Y4C4+Y6C6-Y1C1-Y3C3-Y5C5-Y7C7);

*///在對列做DCTvoididctcol(double*blk)

{

doubletmp[16];

//firststep

tmp[0]=blk[0*8]*C4+blk[2*8]*C2;

tmp[1]=blk[4*8]*C4+blk[6*8]*C6;

tmp[2]=blk[0*8]*C4+blk[2*8]*C6;

tmp[3]=-blk[4*8]*C4-blk[6*8]*C2;

tmp[4]=blk[0*8]*C4-blk[2*8]*C6;

tmp[5]=-blk[4*8]*C4+blk[6*8]*C2;

tmp[6]=blk[0*8]*C4-blk[2*8]*C2;

tmp[7]=blk[4*8]*C4-blk[6*8]*C6;

tmp[8]=blk[1*8]*C7-blk[3*8]*C5;

tmp[9]=blk[5*8]*C3-blk[7*8]*C1;

tmp[10]=blk[1*8]*C5-blk[3*8]*C1;

tmp[11]=blk[5*8]*C7+blk[7*8]*C3;

tmp[12]=blk[1*8]*C3-blk[3*8]*C7;

tmp[13]=-blk[5*8]*C1-blk[7*8]*C5;

tmp[14]=blk[1*8]*C1+blk[3*8]*C3;

tmp[15]=blk[5*8]*C5+blk[7*8]*C7;

//secondstep

tmp[0]=0.5*(tmp[0]+tmp[1]);

tmp[1]=0.5*(tmp[2]+tmp[3]);

tmp[2]=0.5*(tmp[4]+tmp[5]);

tmp[3]=0.5*(tmp[6]+tmp[7]);

tmp[4]=0.5*(tmp[8]+tmp[9]);

tmp[5]=0.5*(tmp[10]+tmp[11]);

tmp[6]=