1、Contents1. DCT變換編碼C語言2. MPEG4中逆DCT變換3. DCT變換探究4. 快速DCT變換DCT變換編碼C語言#include <memory.h>#include <stdio.h>#include <math.h>#include <time.h>#define PI 3.1415926#define CLK_TCK CLOCKS_PER_SECint N;void DCT(double *f,double *F)  int n,m,x; double *dTemp = new double

2、N*N;/中間矩陣 double *coff = new doubleN*N;/變換系數(shù) coff0 = 1/sqrt(N); for( m=1; m<N; m+ )  coffm = sqrt(2)/sqrt(N); memset( dTemp, 0, sizeof(double)*N*N ); memset( F, 0, sizeof(double)*N*N ); /一維變換 for(n=0;n<N;n+)  for(m=0;m<N;m+)  

3、 for(x=0;x<N;x+)    dTempm*N+n += fx*N+n * coffm * cos( (2*x+1) * PI * m/(2*N) ); /第二次一維變換 for(m=0;m<N;m+)c  for(n=0;n<N;n+)   for(x=0;x<N;x+)    Fm*N+n += dTempm*N+x * coffn * cos( (2*x+1) * PI * n/(2*N) );&

4、#160;delete dTemp; delete coff;void iDCT(double *f,double *F) int m,y,x; double *dTemp=new doubleN*N;/中間矩陣 double *coff=new doubleN*N;/變換系數(shù) coff0=1/sqrt(N); for(m=1;m<N;m+)  coffm=sqrt(2)/sqrt(N); memset(dTemp,0,sizeof(double)*N*N); memset(F,0,siz

5、eof(double)*N*N); /一維變換 for(x=0;x<N;x+)  for(y=0;y<N;y+)   for(m=0;m<N;m+)    dTempx*N+y+=Fx*N+m*coffm*cos(2*y+1)*PI*m/(2*N); /第二次一維變換 for(y=0;y<N;y+)  for(x=0;x<N;x+)   for(m=0;m<N;m+) &

6、#160;  Fx*N+y+=dTempm*N+y*coffm*cos(2*x+1)*PI*m/(2*N); delete dTemp; delete coff;int main()   clock_t start,end; start=clock(); int i; long L;  printf("變換維數(shù):"); scanf("%d",&N);  double *f=new doubleN*N

7、;/初始矩陣 double *F=new doubleN*N;/變換后輸出矩陣 memset(F,0,sizeof(double)*N*N);/初始化為0 for(i=0;i<N*N;i+)   printf("f%d%d:",i/N,i%N);  scanf("%lf",&fi);  printf("循環(huán)次數(shù):"); scanf("%d",&L);/輸出初始矩陣 print

8、f("變換前:n"); for(i=1;i<=N*N;i+)   printf("%ft",fi-1);  if(i%N=0)   printf("n");  for(i=0;i<L;i+)  DCT(f,F);/變換/輸出變換后矩陣 printf("變換后:n"); for(i=1;i<=N*N;i+)   print

9、f("%ft",Fi-1);  if(i%N=0)   printf("n");  for(i=0;i<L;i+) iDCT(f,F); /輸出反變換后矩陣 printf("反變換后:n"); for(i=1;i<=N*N;i+)   printf("%ft",fi-1);  if(i%N=0)   printf(&

10、quot;n");  /printf("n"); delete f; delete F; end=clock(); printf("耗時(shí):%fn",(double)(end-start)/CLK_TCK); return 0; =MPEG4中逆DCT變換一旦DCT系數(shù)Fuv被恢復(fù)，那么就可以用逆DCT變換來獲得逆變換值fyx，這些只要被飽和到-256fyx255。對短頭格式，由于不存在隔行模式，因此全部用幀 DCT變換，就是一般的情況。非短頭格式時(shí)，如果使用隔行模式

11、，并且dct_type等于1，此時(shí)使用場DCT變換，場DCT變換的算法同幀DCT變換完全一樣，只是輸出的時(shí)候需要將按場組織的宏塊轉(zhuǎn)換為按幀組織的宏塊。下面簡單介紹一下DCT變換和逆變換的過程。矩陣大小為NxN的二維DCT變換為：u, v, x, y = 0, 1, 2, ¼ N-1其中x, y 是原始域中的空間坐標(biāo)u, v 是變換域中的空間坐標(biāo)逆DCT變換定義為：如果每個(gè)像素為n比特，則DCT變換的輸入為n+1比特，逆DCT變換的輸出為n+1比特。DCT變換后的DCT系數(shù)的長度為n+4比特，動(dòng)態(tài)范圍為-2n+3:+2n+3-1。對我們來說這里的n等于8。NxN的逆DCT變換的實(shí)現(xiàn)必須

12、符合IEEE的關(guān)于8x8的逆DCT變換的實(shí)現(xiàn)的標(biāo)準(zhǔn)，即IEEE Standard Specification for the Implementations of 8 by 8 Inverse Discrete Cosine Transform, Std 1180-1990, December 6, 1990，不過有下列修改：1) IEEE規(guī)范中的3.2小節(jié)的item(1)中最后一句話被替換為：<<Data sets of 1 000 000 (one million) blocks each should be generated for (L=256, H=255),

13、 (L=H=5) and (L=384, H=383). >>2) IEEE規(guī)范中的3.3小節(jié)的text被替換為：<<For any pixel location, the peak error shall not exceed 2 in magnitude. There is no other accuracy requirement for this test.>>3) Let F be the set of 4096 blocks Biyx (i=0.4095) defined as follows :a) Bi00 = i - 2048b) Bi77

14、 = 1 if Bi00 is even, Bi77 = 0 if Bi00 is oddc) All other coefficients Biyx other than Bi00 and Bi77 are equal to 0For each block Biyx that belongs to set F defined above, an IDCT that claims to be compliant shall output a block fyx that as a peak error of 1 or less compared to the reference saturat

15、ed mathematical integer-number IDCT fíí(x,y). In other words, | fyx - fíí(x,y)| shall be <= 1 for all x and y.NOTE 1 lause  2.3 Std 1180-1990 “Considerations of Specifying IDCT Mismatch Errors” requires the specification of periodic intra-picture coding in order to contro

16、l the accumulation of mismatch errors. Every macroblock is required to be refreshed before it is coded 132 times as predictive macroblocks. Macroblocks in B-pictures (and skipped macroblocks in P-pictures) are excluded from the counting because they do not lead to the accumulation of mismatch errors

17、. This requirement is the same as indicated in 1180-1990 for visual telephony according to ITU-T Recommendation H.261.NOTE 2 Whilst the IEEE IDCT standard mentioned above is a necessary condition for the satisfactory implementation of the IDCT function it should be understood that this is not suffic

18、ient. In particular attention is drawn to the following sentence from subclause 5.4: “Where arithmetic precision is not specified, such as the calculation of the IDCT, the precision shall be sufficient so that significant errors do not occur in the final integer values.”逆DCT變換的過程這里不再詳述，需要實(shí)現(xiàn)這個(gè)的可以去參考這

19、個(gè)標(biāo)準(zhǔn)。在實(shí)際應(yīng)用中一般通過兩次1-D IDCT變換來完成2-D IDCT變換，這種方法通常被稱為行列法。一般來說，后者在結(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變換。 圖

20、 48 折疊結(jié)構(gòu)的二維IDCT單元一行數(shù)據(jù) ( 一行有8個(gè)像素?cái)?shù)據(jù) ) 在該單元中的處理流程是：>>>>>>>。DCT變換探究 1 前言 此文適合于那些對DCT或?qū)aar小波的Mallat算法有一定了解的人。 由于我還是高一新丁，文學(xué)底子很薄弱，對于一些技術(shù)方面的知識，我是有口說不出，無法用文字表達(dá)出來，因此這里提供的知識只是我 所知道的1/4左右，還有3/4我不知該如何表達(dá)，特別是第三節(jié)“ 深入研究DCT”，我個(gè)人認(rèn)為簡直是淺入！ 如果你只是菜鳥，不但想看懂此文，而且還要看懂其他的類似文章，那么我教你一個(gè)最快的學(xué)習(xí)方法： 設(shè) X=10,20 分解的

21、方法：低頻=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 只要搞清楚低頻和高頻是怎么來的和如何合并的即可。 2 DCT簡介 DCT全名為Discrete Cosine Transform，中文名為離散余弦變換。在眾人皆知的JPEG編碼中，就是使用了DCT來壓縮圖像的。為什么DCT可以壓縮圖像？我想這個(gè)問題有很多人都想知道，但其實(shí)這是錯(cuò)誤的說法！因?yàn)镈CT在圖像壓縮中僅僅起到扶助的作用，給它n個(gè)數(shù)據(jù)，經(jīng)變換后仍然會(huì)得出n個(gè)數(shù)據(jù)，DCT

22、只不過消除了這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)，就是不

23、能無損變換，因?yàn)镈CT系數(shù)都是一些無理數(shù)，目前為止，依然無法解決。3 深入研究首先讓我們來看看AAN算法的第一階級變換代碼： For I = 0 To 3 J = 7 - I Y(I) = X(I) + X(J) Y(J) = X(I) - X(J) Next I 設(shè)X=10,20,30,40,50,60,70,80 那么Y=90,90,90,90,-10,-30,-50,-70 可以看出，這一階級的低頻部分(相加得出的數(shù)據(jù))全部相等，而高頻部分則呈線性或者是有規(guī)律的。DCT 之所以能以較少的數(shù)據(jù)大致還原圖像，就是因?yàn)橥ㄟ^預(yù)測高頻部分而達(dá)到的。那么為何高頻部分可以預(yù)測呢？請仔細(xì)看上面的代碼，可

24、以 看出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ù)。第二階級變換則是在第一階級變換的基礎(chǔ)上 再次分解出低、高頻，和預(yù)測高頻，得出修正值。第三階級。最后，再將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次之，以此

25、類推?？梢园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)陣的），對于

26、一幅圖像，必須將其分割成無數(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原因暫不討論） 5 DCT的改進(jìn)下面的代碼是我對DCT變換的改進(jìn)，它具有以下特性 l 無損變換 l 計(jì)算量少 l

27、 原位計(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 '-第一階級 For I = 0 To 3 J = 7 - I Y(J) = X(J) - X(I) Y(I) = X(I) + Fix(Y(J) / 2) Next I '-第二階級 For H = 0 To 4 Step 4 For I = 0 To 1 J = 3 - I X(J + H) = Y(J + H) - Y(I + H) X(

28、I + H) = Y(I + H) + Fix(X(J + H) / 2) Next I Next H '-第三階級 For I = 0 To 6 Step 2 Y(I + 1) = X(I + 1) - X(I) Y(I) = X(I) + Fix(Y(I + 1) / 2) Next I '-預(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ì)文平，郭寶龍

29、 Visual C+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)作為下一級節(jié)點(diǎn)，依次進(jìn)行多次，最后得到變換結(jié)果。一、DCT部分例子：Define cos(n*pi/16) CnF(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)

30、*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             

31、60;+x(3)-x(4)*C7從上面的式子可以看到07,16,25,34可以作為第一次運(yùn)算的相加節(jié)點(diǎn)，將所有節(jié)點(diǎn)的表達(dá)式列出后，可發(fā)現(xiàn)一個(gè)規(guī)律，得到一蝶形圖，按之編程，如下：#define  C1 0.9808#define  C2 0.9239#define  C3 0.8315#define  C4 0.7071#define  C5 0.5556#define  C6 0.3827#define  C7 0.1951/先做行DCTvoid fdctrow(double *blk) double S07,S1

32、6,S25,S34,S0734,S1625;double D07,D16,D25,D34,D0734,D1625;S07=blk0+blk7;S16=blk1+blk6;S25=blk2+blk5;S34=blk3+blk4;S0734=S07+S34;S1625=S16+S25;D07=blk0-blk7; D16=blk1-blk6;D25=blk2-blk5;D34=blk3-blk4;D0734=S07-S34;D1625=S16-S25;blk0=0.5*(C4*(S0734+S1625);blk1=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);blk2=0.5

33、*(C2*D0734+C6*D1625);blk3=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);blk4=0.5*(C4*(S0734-S1625);blk5=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);blk6=0.5*(C6*D0734-C2*D1625);blk7=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);/再做列DCTvoid fdctcol(double *blk)double S07,S16,S25,S34,S0734,S1625;double D07,D16,D25,D34,D0734,D1625;S07=

34、blk0*8+blk7*8;S16=blk1*8+blk6*8;S25=blk2*8+blk5*8;S34=blk3*8+blk4*8;S0734=S07+S34;S1625=S16+S25;D07=blk0*8-blk7*8; D16=blk1*8-blk6*8;D25=blk2*8-blk5*8;D34=blk3*8-blk4*8;D0734=S07-S34;D1625=S16-S25;blk0*8=0.5*(C4*(S0734+S1625);blk1*8=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);blk2*8=0.5*(C2*D0734+C6*D1625);bl

35、k3*8=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);blk4*8=0.5*(C4*(S0734-S1625);blk5*8=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);blk6*8=0.5*(C6*D0734-C2*D1625);blk7*8=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);void fdct(double *block)  int i;  for (i=0; i<8; i+)    fdctrow(block+8*i);  for (i=0;

36、 i<8; i+)    fdctcol(block+i);二、IDCT部分圖片來源：G. G. Pechanek, C. W. Kurak, C. J. Glossner, C. H. L. Moller, and S. J. Walsh IBM Microelectronics Division, Research Triangle Park, N.C. "M.F.A.S.T.: A HIGHLY PARALLEL SINGLE CHIP DSP WITH A 2D IDCT EXAMPLE"圖中未給出系數(shù)，需自行算出，仿上述DCT方

37、法直接展開表達(dá)式搜尋規(guī)律即可。編程如下：#define  C1 0.9808#define  C2 0.9239#define  C3 0.8315#define  C4 0.7071#define  C5 0.5556#define  C6 0.3827#define  C7 0.1951/對行做DCTvoid idctrow(double *blk) double tmp16; /first step tmp0=blk0*C4+blk2*C2; tmp1=blk4*C4+blk6*

38、C6; tmp2=blk0*C4+blk2*C6; tmp3=-blk4*C4-blk6*C2; tmp4=blk0*C4-blk2*C6; tmp5=-blk4*C4+blk6*C2; tmp6=blk0*C4-blk2*C2; tmp7=blk4*C4-blk6*C6; tmp8=blk1*C7-blk3*C5; tmp9=blk5*C3-blk7*C1; tmp10=blk1*C5-blk3*C1; tmp11=blk5*C7+blk7*C3; tmp12=blk1*C3-blk

39、3*C7; tmp13=-blk5*C1-blk7*C5; tmp14=blk1*C1+blk3*C3; tmp15=blk5*C5+blk7*C7; /second step tmp0=0.5*(tmp0+tmp1); tmp1=0.5*(tmp2+tmp3); tmp2=0.5*(tmp4+tmp5); tmp3=0.5*(tmp6+tmp7); tmp4=0.5*(tmp8+tmp9); tmp5=0.5*(tmp10+tmp11); tmp6=0.5*(tmp12+tmp13)

40、; tmp7=0.5*(tmp14+tmp15); /third step blk0=tmp0+tmp7; blk1=tmp1+tmp6; blk2=tmp2+tmp5; blk3=tmp3+tmp4; blk4=tmp3-tmp4; blk5=tmp2-tmp5; blk6=tmp1-tmp6; blk7=tmp0-tmp7; /* blk0=0.5*(Y0C4+Y2C2+Y4C4+Y6C6+Y1C1+Y3C3+Y5C5+Y7C7); blk1=0.5*(Y0C4

41、+Y2C6-Y4C4-Y6C2+Y1C3-Y3C7-Y5C1-Y7C5); blk2=0.5*(Y0C4-Y2C6-Y4C4+Y6C2+Y1C5-Y3C1+Y5C7+Y7C3); blk3=0.5*(Y0C4-Y2C2+Y4C4-Y6C6+Y1C7-Y3C5+Y5C3-Y7C1); blk4=0.5*(Y0C4-Y2C2+Y4C4-Y6C6-Y1C7+Y3C5-Y5C3+Y7C1); blk5=0.5*(Y0C4-Y2C6-Y4C4+Y6C2-Y1C5+Y3C1-Y5C7-Y7C3); blk6=0.5*(Y0C4+Y2C6-Y4C4-Y6C2

42、-Y1C3+Y3C7+Y5C1+Y7C5); blk7=0.5*(Y0C4+Y2C2+Y4C4+Y6C6-Y1C1-Y3C3-Y5C5-Y7C7); */在對列做DCTvoid idctcol(double *blk) double tmp16; /first step tmp0=blk0*8*C4+blk2*8*C2; tmp1=blk4*8*C4+blk6*8*C6; tmp2=blk0*8*C4+blk2*8*C6; tmp3=-blk4*8*C4-blk6*8*C2; tmp4=blk0*8*C4-blk2*8*C6; tmp5=-blk4*8*C4+blk6*8*C2; tmp6=blk0*8*C4-blk2*8*C2;