探討JPEG：IDCT

 

轉載 sohu blog:編程鋪路石，作者：shepherd

http://zqw100.blog.sohu.com/6110148.html

 

 

jpeg解碼的核心問題是IDCT（逆離散餘弦變換），這也是整個JPEG解碼過程中耗時最長的步驟。本文將討論一個8×8矩陣的IDCT簡化的問題。
首先看一下IDCT的公式：
   
          7   7                                     2*x+1                   
2*y+1
f(x,y) = sum sum alpha(u)*alpha(v)*F(u,v)*cos (------- *u*PI)* cos
(------ *v*PI)
             u=0 v=0                                    
16                       16
其中x,y=0,1...7 
           { 1/sqrt(8) 
(u==0)
alpha(u) = {
           { 1/2       
(u!=0)
上式中F(u,v)代表未經處理資料（即DCT資料）矩陣中的元素(u,v)，
f(x,y)指解出的像素矩陣中的元素(x,y)
那個公式看暈了吧？好等你明白過來我繼續說
……
明白了哈，我接著說。
看到了那個公式，相信很多人的第一反映就是直接展開，用嵌套的FOR迴圈進行運算。對於這種方法，怎麼說呢，也許用主頻為100G的CPU，解碼的速度還可以忍受。不相信，你可以自己寫一下，看那是幾層迴圈嵌套，再計算一下解碼這64個點一共需要運算多少次。現在理解了吧，直接展開的方法是不行滴。
那麼如何簡化呢？首先就是要把這個二維運算，分解成兩個一維運算。換句話說，先對矩陣的每行進行一次一維IDCT，把結果儲存在一個中間數組裡，然後對這個數組的每列進行一次一維IDCT，就得到最終的結果了。如果你數學基礎還可以的話，可以直接變換上面的公式；如果你的數學比較爛，可以用具體的數字進行演算，你會發現在上文那個醜陋的多重嵌套FOR迴圈中，有許多運算是重複進行的。
好休息一會，自己演算一下把
……
現在我們來講解一維IDCT的簡化。首先看公式
       
1    7                         (2*x+1)*i*PI
p(x) = --- * ∑ ｛ c(i) * DCT(i) *
cos ------------ ｝
        2    i=0                           
16
          { 1/sqrt(2)  (i==0)
其中 c(i)={
          { 1         
(i!=0)
現在假設f(t) = cos(PI*t/16)
展開上式可以得到：
P(0) = (1/4) * (
f(0)*1/sqrt(2)*DCT(0) + f(1)*1*DCT(1) + …… + f(7)*1*DCT(7) )
……
P(7) =
(1/4) * ( f(0)*1/sqrt(2)*DCT(0) + f(15)*1*DCT(1) + …… + f(105)*1*DCT(7)
)
暫時總結一下，將t值整理成一個矩陣：
0  1  2  3  4  5  6  7
0  3  6  9  12 15 18
21
0  5  10 15 20 25 30 35
0  7  14 21 28 35 42 49
0  9  18 27 36 45 54
63
0  11 22 33 44 55 66 77
0  13 26 39 52 65 78 91
0  15 30 45 60 75 90
105
很有規律吧？所以要接著簡化啊。
1，f(t) =
cos(t*PI/16)，說明f(t)是餘弦函數，可以利用餘弦函數的周期性減少未知量的個數。
2，由於f(0)*1/sqrt(2) =
1*sqrt(2)/2 = cos(PI/4) = cos(PI*4/16) = f(4)，利用這一特點可以進一步簡化
假設Cn =
cos(n*PI/16) = f(n)，則一維IDCT的結果可以整理如下：
        DCT(0)*  DCT(1)*  DCT(2)* 
DCT(3)*  DCT(4)*  DCT(5)*  DCT(6)*  DCT(7)* 
P(0) =   C4       C1      
C2       C3       C4       C5       C6       C7
P(1) =   C4       C3      
C6      -C7      -C4      -C1      -C2      -C5
P(2) =   C4       C5     
-C6      -C1      -C4       C7       C2       C3
P(3) =   C4       C7     
-C2      -C5       C4       C3      -C6      -C1
P(4) =   C4      -C7     
-C2       C5       C4      -C3      -C6       C1
P(5) =   C4      -C5     
-C6       C1      -C4      -C7       C6      -C3
P(6) =   C4      -C3      
C6       C7      -C4       C1      -C2       C5
P(7) =   C4      -C1      
C2      -C3       C4      -C5       C6     
-C7
這裡C1-C7都已經是常數了，可以事先求出。
看花眼了吧？好休息一下接著說
……
仔細觀察上面的結果，可以找到規律：
偶數列的上四行和下四行是對稱的；奇數行的上四行和下四行是反對稱的。
假設原始數組是x[]，轉化後的數組是y[],那麼
a0
= x[0]*C4 + x[2]*C2 + x[4]*C4 + x[6]*C6
a1 = x[0]*C4 + x[2]*C6 - x[4]*C4 -
x[6]*C2
a2 = x[0]*C4 - x[2]*C6 - x[4]*C4 + x[6]*C2
a3 = x[0]*C4 - x[2]*C2
+ x[4]*C4 - x[6]*C6

b0 = x[1]*C1 + x[3]*C3 + x[5]*C5 + x[7]*C7
b1 = x[1]*C3 - x[3]*C7 -
x[5]*C1 - x[7]*C5
b2 = x[1]*C5 - x[3]*C1 + x[5]*C7 + x[7]*C3
b3 = x[1]*C7
- x[3]*C5 + x[5]*C3 - x[7]*C1

y[0] = a0 + b0
y[7] = a0 - b0
y[1] = a1 + b1
y[6] = a1 - b1
y[2]
= a2 + b2
y[5] = a2 - b2
y[3] = a3 + b3
y[4] = a3 -
b3
以上是別人文章給出的代碼，實際上還可以進一步化簡a部分
t0 = x[0] * C4
t1 = x[2] * C2
t2 =
x[2] * C6
t3 = x[4] * C4
t4 = x[6] * C6
t5 = x[6] * C2

a0 = t0 + t1 + t3 + t4
a1 = t0 + t2 - t3 - t5
a2 = t0 - t2 - t3 +
t5
a3 = t0 - t1 + t3 - t4
後面的B部分和Y部分不用再變了。

至此化簡工作結束。這樣一個8*8的矩陣，只需進行兩次類似的四則運算就可以了。
我統計了一下，每次一維IDCT是22次乘法，32次加減法，而且經過處理後都是整數運算。
相信你都暈乎乎了吧，好在都說完了，嘿嘿。